Rwanda Baseline Soil Health Study Report

1 Objective

Summary of the demographics and soil characteristics of the Rwanda long term soil health study.

1.1 Study methodology

Add in details and links on study methodology here.

library(knitr)

package ‘knitr’ was built under R version 3.3.2

library(ggplot2)

package ‘ggplot2’ was built under R version 3.3.2

library(stringr)

package ‘stringr’ was built under R version 3.3.2

suppressMessages(library(dplyr))
library(sp)

package ‘sp’ was built under R version 3.3.2

suppressMessages(library(rgdal))

package ‘rgdal’ was built under R version 3.3.2

suppressMessages(library(dismo))

package ‘dismo’ was built under R version 3.3.2

suppressMessages(library(stargazer))
library(leaflet)

package ‘leaflet’ was built under R version 3.3.2

library(XML)

package ‘XML’ was built under R version 3.3.2

suppressMessages(library(maptools))

package ‘maptools’ was built under R version 3.3.2

library(automap)
suppressMessages(library(RStata))
suppressMessages(library(fields))

package ‘fields’ was built under R version 3.3.2

library(gstat)

package ‘gstat’ was built under R version 3.3.2

library(htmltools)

package ‘htmltools’ was built under R version 3.3.2

suppressMessages(library(Matching))

package ‘MASS’ was built under R version 3.3.2

library(reshape2)
options("RStata.StataVersion" = 12)
options("RStata.StataPath" = "/Applications/Stata/StataSE.app/Contents/MacOS/stata-se")
#chooseStataBin("/Applications/Stata/StataSE.app/Contents/MacOS/stata-se")

The working directory objects are from the original baseline analysis. I’ve updated this to pull from a stable working directory (5/4/17)

# wd <- "/Users/mlowes/drive/optimized_agronomy/soil/soil_health_study/data/rw baseline"
# dd <- paste(wd, "data", sep = "/")
# od <- paste(wd, "output", sep="/")
# md <- paste(wd, "maps", sep="/")
drive <- "~/drive/r_help/4_output/statistical_test_outputs"
#load data:
# This data is being drawn from the Soil lab repository. It has the baseline data with it
# d <- read.csv(paste(dd, "Rwanda_shs_commcare_soil_data_final.csv", sep="/"),  stringsAsFactors=FALSE)

2 Data Prep

2.1 Combine survey and soil data

Replicate Alex’s and Emmanuel’s merge process using “Identifiers with SSN_final” provided by Emmanuel. I use the RStata package found here

# I'm replicating Alex's do file here.
stata("merge_shs.do")

. * merge raw commcare data with soil database
. 
. * date: 11 july 2016
. 
. clear all
. set more off
. 
. * set directory
. global wd "/Users/mlowes/drive/soil health study/data/rw baseline"
. global dd "/Users/mlowes/drive/soil health study/data/rw baseline/data"
. global troubleshoot "/Users/mlowes/drive/soil health study/data/rw baseline/t
> roubleshoot"
. 
. * insheet data
. insheet using "$dd/raw_rwanda_commcare_shs.csv", clear
file /Users/mlowes/drive/soil health study/data/rw
    baseline/data/raw_rwanda_commcare_shs.csv not found
r(601);
. 
. * clean sample_id variable
. replace sample_id = "" if sample_id == "---"
. replace sample_id = "-99" if sample_id == "99"
. 
. destring sample_id, replace
. 
. * manipulate current "sample_id" to become a proper alpha-numeric unique iden
> tifier
. * this simply involves adding "C" to each numeric code that belongs to a cont
> rol farmer
. 
. * variables of interest:
.       * d_client
.       * sample_id
. 
. ren demographic_infod_client d_client
. replace d_client = "" if d_client == "---"
. destring d_client, replace
. 
. tostring sample_id, gen(sample_id_string)
. 
. replace sample_id_string = sample_id_string + "C" if d_client == 0
. drop sample_id 
. ren sample_id_string sample_id
. 
. ***** !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ***
> ***
. * quantify and note cases where sample_id appears more than twice
.      * 39 codes appear 2x
.      * 1 code appears 3x
.      * 1 code appears 4x
.      * 1 code appears 10x
.      * 1 code appears 19x
. 
. duplicates report sample_id
. duplicates tag sample_id, gen(n_duplicate_id)
. gen d_id_problem = 1 if n_duplicate_id != 0
. replace d_id_problem = 0 if missing(d_id_problem)
. 
. ** need to investigate 5% of observations (114 samples), i.e. d_id_problem = 
> 1**
. ***** !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ***
> ***
. 
. * outsheet problematic observations and examine 1-by-1
. sort sample_id
. outsheet demographicnom_enqueteur district demographic_infocellule_selected d
> emographic_infoumudugudu d_client demographic_infonom_cultivateur sample_id d
> emographic_infosex demographic_infoage_cultivateur demographic_infon_telephon
> e_cult d_id_problem using "$dd/shs_to_check.csv" if d_id_problem == 1, replac
> e comma
. 
. **************** post-investigation correction of incorrect ids *************
> ***
. 
. ren demographic_infonom_cultivateur farmer_name
. ren demographic_infocellule_selected selected_cell
. ren demographicnom_enqueteur enumerator_name
. 
. * drop test/accidentally double-entered observations
. drop if infoformid == "f138897f-6e80-4b0a-9db3-ab2134cd9e51"
. drop if farmer_name == "Mugisha"
. drop if farmer_name == "Muhunde yoramu"
. drop if farmer_name == "Test"
. drop if farmer_name == "Billy"
. drop if farmer_name == "Gakuba michel"
. drop if farmer_name == "Renata"
. drop if farmer_name == "Bub"
. drop if farmer_name == "M"
. drop if farmer_name == "Anita"
. drop if farmer_name == "Jean"
. drop if farmer_name == "Nyirazaninka lea"
. drop if farmer_name == "Mugiraneza Pacifiue"
. drop if farmer_name == "Nsengiyaremye Innocent"
. 
. * correct incorrectly recorded ids
. replace sample_id = "880C" if farmer_name == "Nyiransanzineza Marie Rose"
. replace sample_id = "811C" if farmer_name == "Ngirabatware  Daniel"
. replace sample_id = "646C" if farmer_name == "Nyirasabwa anonciata"
. replace sample_id = "624C" if farmer_name == "Nyirinkindi j Bosco"
. replace sample_id = "560C" if farmer_name == "Mubirigi Bertin"
. replace d_client = 0 if farmer_name == "Mubirigi Bertin"
. replace sample_id = "375C" if farmer_name == "Hitimana ezecquiere"
. replace sample_id = "322C" if farmer_name == "Nyirantabire arrivera"
. replace sample_id = "2566C" if farmer_name == "MUSABYIMANA Colletta"
. replace sample_id = "2399C" if farmer_name == "Nyirabajyambere anastisia"
. replace d_client = 0 if farmer_name == "Nyirabajyambere anastisia"
. replace sample_id = "2280C" if farmer_name == "Zihabake Simon"
. replace sample_id = "2202C" if farmer_name == "mugiraneza pacifique"
. replace sample_id = "1780C" if farmer_name == "Hakuzimana Theophile"
. replace sample_id = "1757C" if farmer_name == "Bavugirije Feleciane"
. replace sample_id = "1626C" if farmer_name == "Nyabivumu"
. replace sample_id = "1575C" if farmer_name == "BAHIMBANDE Beriyana."
. replace d_client = 0 if farmer_name == "BAHIMBANDE Beriyana."
. replace sample_id = "1103C" if farmer_name == "Barirwanda Elyse"
. replace sample_id = "1037C" if farmer_name == "Mukashema Faraziya"
. replace d_client = 0 if farmer_name == "Mukashema Faraziya" 
. replace sample_id = "3069" if farmer_name == "Ndababonye Silas"
. replace sample_id = "2611" if farmer_name == "Sebazungu modeste"
. replace d_client = 1 if farmer_name == "Sebazungu modeste"
. replace sample_id = "2415" if farmer_name == "Nkaka joseph"
. replace sample_id = "2412" if farmer_name == "Ndagijimana alexis"
. replace sample_id = "2410" if farmer_name == "Sebuhoro elie"
. replace sample_id = "2406" if farmer_name == "Mugisha ruth"
. replace sample_id = "2405" if farmer_name == "Mboneza hasan"
. replace sample_id = "2404" if farmer_name == "Ntibitonderwa veriane"
. replace sample_id = "2399" if farmer_name == "Nyiraburakeye feresita"
. replace sample_id = "2395" if farmer_name == "Uwimana danier"
. replace sample_id = "2394" if farmer_name == "Bizinde silas"
. replace sample_id = "2390" if farmer_name == "Nyirahabimvana keresesiya"
. replace sample_id = "2388" if farmer_name == "Nyirasekerabanzi tharcilla"
. replace sample_id = "2386" if farmer_name == "Ndererimana emmanuel"
. replace sample_id = "2291" if farmer_name == "Ndacyayisenga Feresiyani"
. replace sample_id = "2184" if farmer_name == "Muhawenimana Beretirida"
. replace sample_id = "2145" if farmer_name == "Fapfakwita celestine"
. replace sample_id = "2141" if farmer_name == "Nsangiranabo krizoroji"
. replace sample_id = "2140" if farmer_name == "Nsabyumuremyi anakereti"
. replace sample_id = "2136" if farmer_name == "Nyirantibitonderwa savoroniya"
. replace sample_id = "2131" if farmer_name == "Hagenimana Ewufroni"
. replace sample_id = "2066" if farmer_name == "Hitayezu vicent"
. replace d_client = 1 if farmer_name == "Hitayezu vicent"
. replace sample_id = "2044" if farmer_name == "Niyodushima janette"
. replace d_client = 1 if farmer_name == "Niyodushima janette"
. replace sample_id = "2024" if farmer_name == "Mukabatsinda veren A"
. replace sample_id = "2020" if farmer_name == "Sibomana innocent"
. replace sample_id = "2009" if farmer_name == "Nsabimana inyasi"
. replace sample_id = "2002" if farmer_name == "Mateso elizabet"
. replace sample_id = "1889" if farmer_name == "Niyonsaba Daniel"
. replace sample_id = "1864" if farmer_name == "Mukankubana Souzanne"
. replace sample_id = "1780" if farmer_name == "Nyirantihabose Annonciatha"
. replace sample_id = "1771" if farmer_name == "Bizimungu Damascne"
. replace sample_id = "1675" if farmer_name == "Bikorimana Isaac"
. replace sample_id = "1103" if farmer_name == "Munyabarata japhette"
. replace sample_id = "954" if farmer_name == "Hahirwabake silver"
. replace sample_id = "541" if farmer_name == "Rugemintwaza Vianney"
. replace sample_id = "454" if farmer_name == "Niyonizeye fororida"
. replace sample_id = "407" if farmer_name == "Munyaneza jean pierre"
. replace sample_id = "375" if farmer_name == "Munyarubuga nanias"
. replace sample_id = "322" if farmer_name == "Gacandaga zackarie"
. replace sample_id = "262" if farmer_name == "MUHIRWA J Damascene"
. replace sample_id = "32" if farmer_name == "Nsanzemuhire"
. 
. *** re-run duplicates test
. drop d_id_problem
. drop n_duplicate_id
. 
. duplicates report sample_id
. duplicates tag sample_id, gen(n_duplicate_id)
. gen d_id_problem = 1 if n_duplicate_id != 0
. replace d_id_problem = 0 if missing(d_id_problem)
. 
. * second try: outsheet problematic observations and examine 1-by-1
. sort sample_id
. outsheet using "$dd/shs_to_check_2.csv" if d_id_problem == 1, replace comma
. 
. * second try: correct incorrect ids
. 
. replace sample_id = "2202" if farmer_name == "Mukamurara scholastique"
. replace d_client = 1 if farmer_name == "Mukamurara scholastique"
. replace sample_id = "2195" if farmer_name == "Ndacyayisenga Feresiyani"
. replace sample_id = "824" if farmer_name == "Ngirabatware  Daniel"
. 
. replace sample_id = "" if rownumber == 238
. replace sample_id = "" if rownumber == 848
. replace sample_id = "" if rownumber == 1257
. replace sample_id = "" if rownumber == 1732
. 
. *** re-run duplicates test
. drop d_id_problem
. drop n_duplicate_id
. 
. duplicates report sample_id
. duplicates tag sample_id, gen(n_duplicate_id)
. gen d_id_problem = 1 if n_duplicate_id != 0
. replace d_id_problem = 0 if missing(d_id_problem)
. 
. * third try: outsheet problematic observations
. sort sample_id
. outsheet using "$dd/shs_to_check_3.csv" if d_id_problem == 1, replace comma
. 
. *** 4 sets of duplicates remain (2 observations)
. drop if d_id_problem == 1
. save "$dd/rwanda_commcare_shs_clean.dta", replace
. outsheet using "$dd/rwanda_commcare_shs_clean.csv", comma replace
. 
. * merge datasets
. insheet using "$dd/mne_rwanda_database_shs.csv", clear
. 
. save "$dd/mne_rwanda_database_shs.dta", replace
. 
. use "$dd/rwanda_commcare_shs_clean.dta", clear
. 
. merge 1:1 sample_id using "$dd/mne_rwanda_database_shs.dta"
. 
.    * 93 observations unmatched... :(
. 
. * investigate unmatched observations
. sort sample_id
. outsheet using "$troubleshoot/unmatched_commcare.csv" if _merge == 1, comma r
> eplace
. outsheet using "$troubleshoot/unmatched_inventory.csv" if _merge == 2, comma 
> replace
. 
. * make first round of corrections in master data
. use "$dd/rwanda_commcare_shs_clean.dta", clear
. 
. replace sample_id = "2711C" if farmer_name == "Kamuhanda Elie"
. replace sample_id = "260C" if farmer_name == "HATEGEKIMANA Mathias"
. replace sample_id = "245" if farmer_name == "Kayiranga Michel"
. replace sample_id = "1543" if farmer_name == "Ndagijimana Eliezel"
. replace sample_id = "1543C" if farmer_name == "Mukantanganzwa odeta"
. replace sample_id = "1632" if farmer_name == "Nyiranziguye Cecile"
. replace sample_id = "2782" if farmer_name == "Mukamungu leocadie"
. replace sample_id = "421C" if farmer_name == "Niyonizeye Francine"
. replace sample_id = "828C" if farmer_name == "Mukamazimpaka Terese"
. replace sample_id = "364" if farmer_name == "Mukamuyango verina"
. replace sample_id = "364C" if farmer_name == "Mujawimana jeanne"
. replace sample_id = "368" if farmer_name == "Nyiramanzi jean damascene"
. replace sample_id = "368C" if farmer_name == "Karwiyegura rachel"
. replace sample_id = "391" if farmer_name == "Mugemane augistin"
. replace sample_id = "391C" if farmer_name == "Ngirishuti sumayire"
. replace sample_id = "395" if farmer_name == "Niyomugabo amiel"
. replace sample_id = "395C" if farmer_name == "Ntawiniga augustin"
. replace sample_id = "396" if farmer_name == "Mpayimana philippe"
. replace sample_id = "329" if farmer_name == "Mucyezangango emmanuel"
. replace sample_id = "411" if farmer_name == "Nsabimana mathias"
. replace sample_id = "411C" if farmer_name == "Singirankabo"
. replace sample_id = "414" if farmer_name == "Rutaharama eldephonse"
. replace sample_id = "329C" if farmer_name == "Mukamana josephine"
. replace sample_id = "657" if farmer_name == "Nyandwi vincent"
. replace sample_id = "569C" if farmer_name == "Uwimana Bercilla"
. replace sample_id = "2051C" if farmer_name == "Ntambabazi Anastase"
. replace sample_id = "2249" if farmer_name == "Uwineza valentine"
. replace sample_id = "2360C" if farmer_name == "Bamporineza j deDieu"
. replace sample_id = "2001" if farmer_name == "Nyiramakuba thanene"
. replace sample_id = "2897C" if farmer_name == "Mukagatanazi Marianne"
. replace sample_id = "2965C" if farmer_name == "Twizeyimana Theoneste"
. 
. drop if farmer_name == "Yjk"
. drop if farmer_name == "Uwambajimana Agnes"
. drop if farmer_name == "Tr"
. 
. replace sample_id = "2415C" if farmer_name == "Nyiramahirwe Frolance"
. replace sample_id = "851" if farmer_name == "Bigirimana  Amon"
. 
. 
. * duplicates check
. drop d_id_problem
. drop n_duplicate_id
. 
. duplicates report sample_id
. duplicates tag sample_id, gen(n_duplicate_id)
. gen d_id_problem = 1 if n_duplicate_id != 0
. replace d_id_problem = 0 if missing(d_id_problem)
. 
. * second try: merge
. merge 1:1 sample_id using "$dd/mne_rwanda_database_shs.dta"
. 
. * deal with 30 unmatched cases
. sort sample_id
. outsheet using "$troubleshoot/unmatched_commcare_2.csv" if _merge == 1, comma
>  replace
. outsheet using "$troubleshoot/unmatched_inventory_2.csv" if _merge == 2, comm
> a replace
. 
. * 5 submitted surveys with no corresponding sample
. *drop if _merge == 1
. 
. * 19 samples with no corresponding survey
. *drop if _merge == 2
. 
. * outsheet merged database
. outsheet using "$dd/rwanda_shs_baseline_data.csv", comma replace
. * THIS IS THE CLEAN VERSION! ONLY USE THIS! DON'T USE OTHER THINGS!
. 

Import the results of Alex’s do file.

d <- read.csv("data/rwanda_shs_baseline_data.csv", stringsAsFactors = F)
Identifiers <- read.csv("data/Identifiers with SSN_final.csv", stringsAsFactors = F)
wetChem <- read.csv(paste("/Users/mlowes/drive/jupyter/Rwanda Acidity", "Original chem_rwanda_shs.csv", sep = "/"))

d <- left_join(d, Identifiers, by="sample_id")
# now import the soil results and merge with survey data
soilRF <- read.csv("data/rwshs_rfresults.csv", stringsAsFactors = F)
names(soilRF)[names(soilRF)=="X"] <- "SSN"
d <- left_join(d, soilRF, by="SSN")

2.2 Cleaning baseline variables

Now let’s start cleaning the demographic variables

d[d=="---"] <- NA
names(d) <- gsub("text_final_questions", "", names(d))
names(d) <- gsub("intro_champ_echantillon", "", names(d))
names(d) <- gsub("demographic_info", "", names(d))
names(d) <- gsub("other_inputs_", "", names(d))
names(d) <- gsub("crop1_15b_inputs", "", names(d))
names(d) <- gsub("crop2_15b_inputs", "", names(d))
names(d) <- gsub("^15b", "", names(d))
names(d) <- gsub("historical_intro", "", names(d))
names(d)[names(d)=="field_dim"] <- "field_dim1"
names(d)[names(d)=="v51"] <- "field_dim2"

Address unusual field sizes

ggplot(d, aes(x=field_dim1, y=field_dim2)) + geom_point()

It seems really unlikely that fields are 200 meters long while only being 20 meters wide. I don’t know how to check this though.

# clean field dimensions here - winsor the values to something reasonable.
d[(d$field_dim1>=100 | d$field_dim2>=100) & !is.na(d$field_dim1), c("field_dim1", "field_dim2")]

Take care of demographic data formatting issues

# deal with names and drop unnecessary variables
d <- d %>% 
    dplyr::select(-c(rownumber, infoformid, introductiond_accept, photo,
        infocompleted_time, 
        enumerator_name, contains("phone"), farmer_name, farmersurname, farmername,
        d_respondent, additionalsamplepackedandsenttol, additionalsamplerequestedfromlab,
        datedryingcompleteifnecessary, driedindistrictifnecessary, senttohqyo,
        collectedindistrictyo, excessstoredathq_, receivedathq_,dateofinitialdryingifnecessary,
        samplecollectedinfieldyo, field_des, samplewetordry)) %>%
    rename(
    female = sex,
    age = age_cultivateur,
    own = d_own,
    client = d_client) %>%
    mutate(
    female= ifelse(female=="gore", 1,0),
    field.size = field_dim1*field_dim2
    )
d$total.seasons <- apply(d[, grep("d_season_list", names(d))], 1, function(x) {
    sum(x, na.rm=T)})

Clean agronomic practice variables

agVars <- c("n_season_fert", "n_season_compost", "n_season_lime", "n_season_fallow",
            "n_seasons_leg_1", "n_seasons_leg_2")
summary(d[,agVars])

 n_season_fert    n_season_compost n_season_lime     n_season_fallow   n_seasons_leg_1  n_seasons_leg_2   
 Min.   : 0.000   Min.   : 0.000   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000   Min.   : -88.000  
 1st Qu.: 0.000   1st Qu.: 2.000   1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:   0.000  
 Median : 1.000   Median : 5.000   Median : 0.0000   Median : 0.0000   Median : 1.000   Median :   2.000  
 Mean   : 2.041   Mean   : 5.651   Mean   : 0.1819   Mean   : 0.6355   Mean   : 2.171   Mean   :   4.368  
 3rd Qu.: 3.000   3rd Qu.:10.000   3rd Qu.: 0.0000   3rd Qu.: 0.0000   3rd Qu.: 4.000   3rd Qu.:   5.000  
 Max.   :10.000   Max.   :10.000   Max.   :10.0000   Max.   :10.0000   Max.   :10.000   Max.   :3333.000  
 NA's   :19       NA's   :19       NA's   :19        NA's   :19        NA's   :19       NA's   :19        

Sort out the legumes as a second crop

table(d$n_seasons_leg_2, useNA = 'ifany')

 -88    0    1    2    3    4    5    6    7    8    9   10   15   22   50   53   83   88  101  102 3333 <NA> 
   1  956  221  258  174  130  301   61   34   66   57  199    1    1    1    3    1    1    1    1    1   19 

d$n_seasons_leg_2 <- ifelse(d$n_seasons_leg_2 <0 | d$n_seasons_leg_2>10, NA, d$n_seasons_leg_2)

2.2.1 Check out client variable - why are there missing values?

table(d$client, useNA = 'ifany')

   0    1 <NA> 
1233 1236   19 

d[is.na(d$client), c("sample_id")]

 [1] "1189C" "1207C" "1295"  "1295C" "1298C" "1300"  "1300C" "1366C" "1476C" "1485C" "1554"  "1573"  "1614"  "2042"  "2371C" "2373"  "2375" 
[18] "2382"  "2442" 

# replace client based on whether there is a C in the client variable.
d$client.check <- ifelse(grepl("C", d$sample_id)==T, 0, 1)
table(d$client, d$client.check)

   
       0    1
  0 1233    0
  1    1 1235

It looks like most farmers were recorded correctly except for one farmer who was coded as Tubura farmer but their sample_id indicate their a control. Let’s take a look:

d[d$client==1 & d$client.check==0 & !is.na(d$d_gps), c("sample_id", "tuburacontroltc")]

# they should be a control
d[d$client==1 & d$client.check==0 & !is.na(d$d_gps), "client"] <- 0
# remove farmers for which we have soil data but no survey data (using)
d <- d[-which(grepl("using", d$X_merge)),]
table(d$client, d$client.check, useNA = 'ifany')

   
       0    1
  0 1234    0
  1    0 1235

# update client to equal client check
d$client <- d$client.check

Fix some more variable names:

names(d)[names(d)=="field_kg_fert1_1"] <- "field_kg_fert1_15b"
names(d)[names(d)=="field_kg_fert2_1"] <- "field_kg_fert2_15b"
names(d)[names(d)=="field_kg_compost"] <- "field_kg_compost_15b"

Recode variables to numeric:

# recode to numeric
varlist <- c("client", "own", "crop1_15b_seedkg", "crop1_15b_yield", "crop1_15b_yield_",
    "crop2_15b_seedkg", "crop2_15b_yield", "crop2_15b_yield_", "field_kg_fert1_15b",
    "field_kg_fert2_15b", "field_kg_compost_15b", "d_lime_15b", "kg_lime_15b")
# check that there aren't values hidden in the character variables
#apply(d[,varlist], 2, function(x){table(x, useNA='ifany')})
# recode characters to numerics
d[, varlist] <- sapply(d[,varlist], as.numeric)

2.2.2 Sort out kgs of lime applied

table(d$kg_lime_15b, useNA = 'ifany')

 -88    1    3    4    5   10   13   15   20   25   30   35   50   60   88  100  150  200 <NA> 
   1    2    2    2    6    5    1    1    2    5    1    1    8    1    2    2    4    1 2422 

d$kg_lime_15b <- ifelse(abs(d$kg_lime_15b)==88, NA, d$kg_lime_15b)

# divide out GPS coordinates
# http://rfunction.com/archives/1499
# replace the blank gps_pic_guide with info
d <- cbind(d, str_split_fixed(d$gps_pic_guid, " ", n=4))
names(d)[names(d)=="1" |names(d)== "2" | names(d)== "3" | names(d)== "4"] <- c("lat", "lon", "alt", "precision")
d[,c("lat", "lon", "alt", "precision")] <- sapply(d[,c("lat", "lon", "alt", "precision")],
                                                  function(x){as.numeric(as.character(x))})

2.3 Cleaning soil data

Cleaning of soil data: Come back, check and clean the soil data before outputting to clean data set. Plot each of the soil variables to look for unrealistic values.

dim(d[is.na(d$m3.Ca),])

[1]  26 111

d <- d[-which(is.na(d$m3.Ca)),]
       
summary(d[,c("m3.Ca", "m3.Mg", "pH", "Total.N", "Total.C", "ExAl")])

     m3.Ca            m3.Mg               pH           Total.N           Total.C            ExAl        
 Min.   : 213.8   Min.   :  36.93   Min.   :4.549   Min.   :0.05723   Min.   :0.9107   Min.   :0.02589  
 1st Qu.: 436.5   1st Qu.: 115.94   1st Qu.:5.020   1st Qu.:0.13577   1st Qu.:1.7298   1st Qu.:0.09995  
 Median : 723.1   Median : 184.63   Median :5.535   Median :0.15483   Median :2.1164   Median :0.23170  
 Mean   : 877.1   Mean   : 208.23   Mean   :5.543   Mean   :0.15706   Mean   :2.1147   Mean   :0.57795  
 3rd Qu.:1159.5   3rd Qu.: 270.55   3rd Qu.:6.010   3rd Qu.:0.17910   3rd Qu.:2.3664   3rd Qu.:1.00947  
 Max.   :3724.8   Max.   :1009.88   Max.   :7.131   Max.   :0.24656   Max.   :4.2212   Max.   :2.31905  

2.3.1 Soil Type

names(d)[names(d)=="general_field_infosoil_type"] <- "soil_type"
names(d)[names(d)=="general_field_infotexture_soil"] <- "soil_texture"
d$black_soil <- ifelse(d$soil_type=="black", 1,0)
d$red_soil <- ifelse(d$soil_type=="red", 1,0)
d$sandy_soil <-ifelse(d$soil_type=="sandy", 1,0)

Let’s check out the rows for which we don’t have soil data and drop them as they won’t contribute to the full picture.

2.3.2 Graphs of RW baseline soil variables

soilVars <- c("m3.Ca", "m3.Mg", "pH", "Total.N", "Total.C", "ExAl")

for(i in 1:length(soilVars)){
print(
  ggplot(data=d, aes(x=as.factor(client), y=d[,soilVars[i]])) + 
    geom_boxplot() +
    labs(x="Tubura Farmer", y=soilVars[i], title = paste("RW baseline soil - ", soilVars[i], sep = ""))
  )  
}

2.3.3 Check out soil relationships

There are biologically predictable relationships between soil chemical characteristics. For instance, we expect Ca and Mg to move in the same direction and be positively correlated with pH. If we had Aluminum as an outcome, we’d expect pH to be negatively correlated with soluable aluminum. Let’s look quickly to confirm if those relationships are present:

ggplot(d, aes(x=m3.Ca, y=m3.Mg)) + geom_point() +
    stat_smooth(method="loess") +
    labs(x = "Calcium (m3)", y= "Magnesium (m3)", title="Calcium and Magnesium relationship")

ggplot(d, aes(x=pH, y=m3.Ca)) + geom_point() +
  stat_smooth(method="loess") +
  labs(x = "pH", y="Calcium (m3)", title = "pH and Calcium relationship")

ggplot(d, aes(x=pH, y=m3.Mg)) + geom_point() +
  stat_smooth(method="loess") +
  labs(x = "pH", y="Magnesium (m3)", title = "pH and Magnesium relationship")

ggplot(d, aes(x=pH, y=ExAl)) + geom_point() +
  stat_smooth(method="loess") +
  labs(x = "pH", y="Exchangeable Aluminum", title = "pH and Aluminum relationship")

ggplot(d, aes(x=Total.C, y=Total.N)) + geom_point() + 
  stat_smooth(method="loess") +
  labs(x = "Total Carbon", y="Total Nitrogen", title = "Carbon and Nitrogen relationship")

The soil characteristics are moving in the manner that is consistent with our understanding of soil chemical processes.

Save clean demographic and soil data to external file

write.csv(d, file="data/shs rw baseline.csv")
save(d, file="data/shs rw baseline.Rdata")

2.3.4 Map of baseline observations

Produce a simple map of where our observations are

See here for more on using markerClusterOptions in leaflet.

In the map below, the larger green circles are Tubura farmers and the smaller blue circles are control farmers.

e <- d[!is.na(d$lon),]
ss <- SpatialPointsDataFrame(coords = e[, c("lon", "lat")], data=e)
pal <- colorNumeric(c("navy", "green"), domain=unique(ss$client))
map <- leaflet() %>% addTiles() %>%
  setView(lng=rwanda$longitude, lat=rwanda$latitude, zoom=8) %>%
  addCircleMarkers(lng=ss$lon, lat=ss$lat, 
                   radius= ifelse(ss$client==1, 10,6),
                   color = pal(ss$client),
clusterOptions = markerClusterOptions(disableClusteringAtZoom=13, spiderfyOnMaxZoom=FALSE))
map

3 Summary statistics

3.0.1 Table of final baseline breakdown

3.1 Baseline balance

Let’s see how balanced our farmers are in terms of demographic variables. Tubura farmers were selected based on (list criteria) and control farmers in the same area tha fit the same criteria were also selected. No matching process has been performed to identify the control farmers that most closely resemble the Tubura farmers in the sample. These results are entirely reflecting the balance inherent in the identification process, not any statistical matching of treatment and control.

d$valley <- ifelse(d$general_field_infofield_location=="valley", 1,0)
d$hillside <- ifelse(d$general_field_infofield_location=="hillside", 1,0)
d$hilltop <- ifelse(d$general_field_infofield_location=="hilltop", 1,0)
names(d)[names(d)=="general_field_infograde_hill"] <- "slope"

cor(d[, grep("betail_", names(d))], use='complete.obs')

                       betail_ownedn_inka betail_ownedn_ihene betail_ownedn_inkoko betail_ownedn_ingurube betail_ownedn_intama
betail_ownedn_inka             1.00000000         0.036402152           0.09740274            0.021664199          0.042450521
betail_ownedn_ihene            0.03640215         1.000000000           0.12742469           -0.009667001          0.002625236
betail_ownedn_inkoko           0.09740274         0.127424693           1.00000000            0.049416370          0.030415986
betail_ownedn_ingurube         0.02166420        -0.009667001           0.04941637            1.000000000          0.025487680
betail_ownedn_intama           0.04245052         0.002625236           0.03041599            0.025487680          1.000000000

names(d)[grep("betail_", names(d))] <- c("cows", "goats", "chickens", "pigs", "sheep")
d$own.cows <- ifelse(d$cows>0 & !is.na(d$cows), 1,0)
d$own.goats <- ifelse(d$goats>0 & !is.na(d$goats), 1,0)
d$own.chickens <- ifelse(d$chickens>0 & !is.na(d$chickens), 1,0)
d$own.pigs <- ifelse(d$pigs>0 & !is.na(d$pigs), 1,0)
d$own.sheep <- ifelse(d$sheep>0 & !is.na(d$sheep), 1,0)

3.2 Crop selection in 15B

names(d)[names(d)=="inputuse_15b_priorculture_15b_1"] <-"crop_15b"
d$crop_15b <- tolower(d$crop_15b)
sort(prop.table(table(d$crop_15b, useNA = 'ifany')))

         kor         hwag        sheke       coffee          cer        uburo         teke       insina         shaz          nyo          gan 
0.0004093328 0.0008186656 0.0008186656 0.0012279984 0.0024559967 0.0036839951 0.0040933279 0.0057306590 0.0085959885 0.0135079820 0.0196479738 
   other_veg          ray          yum       ntacyo          gor         soya          jum          big         saka          shy 
0.0212853050 0.0519852640 0.0577159230 0.0589439214 0.0597625870 0.0790012280 0.0920998772 0.1604584527 0.1649611134 0.1927957429 

d$climbing_bean <- ifelse(d$crop_15b=="shy", 1,0)
d$sorghum <- ifelse(d$crop_15b=="saka", 1,0)
d$bush_bean <- ifelse(d$crop_15b=="big", 1,0)
d$maize <- ifelse(d$crop_15b=="gor", 1,0)

out.list <- c("female", "age", "hhsize", "own", "field.size",
              "n_season_fert", "n_season_compost", "n_season_lime", "n_season_fallow", "n_seasons_leg_1", "n_seasons_leg_2", "slope",  "alt", "valley", "hillside", "hilltop", "cows", "goats", "chickens", "pigs", "sheep", "climbing_bean", "sorghum", "bush_bean","maize", soilVars, "own.cows", "own.goats", "own.chickens", "own.pigs", "own.sheep", "black_soil", "red_soil", "sandy_soil")
output <- do.call(rbind, lapply(out.list, function(x) {
  
  out <- t.test(d[,x] ~ d[,"client"], data=d)
  tab <- data.frame(out[[5]][[2]],out[[5]][[1]], out[3])
  tab[,1:2] <- round(tab[,1:2],3)
  names(tab) <- c(names(out[[5]]), "pvalue")
  return(tab)
}))
# use p.adjust with bonferroni correction
output$pvalue <- p.adjust(output$pvalue, method="fdr")
rownames(output) <- out.list
output <- output[order(output$pvalue),]
output$pvalue <- ifelse(output[, 3] < 0.001, "< 0.001", round(output[, 3], 3)) 
colnames(output) <- c("Tubura Client","Non-Tubura", "p-value")  

print(kable(output))

	Tubura Client	Non-Tubura	p-value
n_season_fert	3.263	0.833	< 0.001
female	0.493	0.649	< 0.001
own.cows	0.683	0.539	< 0.001
age	44.021	48.088	< 0.001
hhsize	5.418	4.894	< 0.001
own.chickens	0.342	0.267	< 0.001
n_seasons_leg_2	2.550	3.072	< 0.001
n_season_lime	0.228	0.132	< 0.001
own.pigs	0.355	0.281	< 0.001
pigs	0.525	0.394	0.006
chickens	1.225	0.895	0.007
hillside	0.682	0.733	0.017
own	0.932	0.957	0.025
cows	1.178	0.912	0.025
valley	0.266	0.222	0.028
maize	0.072	0.048	0.028
sorghum	0.146	0.184	0.029
climbing_bean	0.212	0.174	0.039
own.sheep	0.104	0.077	0.044
sheep	0.202	0.142	0.053
n_season_compost	5.818	5.523	0.103
Total.N	0.156	0.158	0.111
n_season_fallow	0.690	0.569	0.111
m3.Mg	204.339	212.129	0.185
m3.Ca	859.917	894.354	0.197
field.size	667.734	601.671	0.232
Total.C	2.104	2.126	0.426
bush_bean	0.154	0.167	0.494
n_seasons_leg_1	2.228	2.135	0.512
hilltop	0.051	0.044	0.512
pH	5.533	5.553	0.512
slope	10.203	10.445	0.592
ExAl	0.584	0.571	0.709
goats	1.076	1.053	0.828
sandy_soil	0.145	0.149	0.828
alt	1570.267	1565.917	0.884
black_soil	0.537	0.541	0.884
red_soil	0.279	0.277	0.93
own.goats	0.485	0.485	0.981

Pre-PSM table 1 summary

table1vars <- c("pH", "Total.C", "Total.N", "m3.Ca", "m3.Mg")
table1 <- output[rownames(output) %in% table1vars, ]
write.csv(table1, file=paste("output","rw table1.csv", sep = "/"), row.names=T)

# write baseline balance table for ES per his table layouts
t4order <- c("age", "female", "hhsize", "own")
table4vars <- paste(t4order, collapse="|")
rw.table4 <- output[grepl(table4vars, rownames(output)),]
rw.table4 <- rw.table4[order(match(rownames(rw.table4), t4order)),]
write.csv(rw.table4, file=paste("output", "pre match balance table4.csv", sep="/"),
          row.names = T)
t5order <- c("field.size", "slope", "alt", "hilltop", "hillside", "valley",
                      "climbing_bean", "sorghum", "bush_bean", "maize")
table5vars <- paste(t5order,collapse="|")
rw.table5 <- output[grepl(table5vars, rownames(output)),]
rw.table5 <- rw.table5[order(match(rownames(rw.table5), t5order)),]
write.csv(rw.table5, file=paste("output", "pre match balance table5.csv", sep = "/"),
          row.names = T)
t6order <- c("own.cows", "cows", "own.pigs", "pigs", "own.sheep", "sheep", "own.goats", "goats","own.chickens", "chickens")
table6vars <- paste(t6order, collapse = "|")
rw.table6 <- output[grep(table6vars, rownames(output)), ]
rw.table6 <- rw.table6[order(match(rownames(rw.table6), t6order)),]
write.csv(rw.table6, file=paste("output", "pre match balance table6.csv", sep = "/"), 
          row.names=T)

3.3 Overall balance interpretation

Demographic variables We are not well balanced along the main demographic variables we collected, sex, age and HH size. For the purposes of inference we can test some matching algorithms to improve the match between Tubura and control farmers.

Agricultural practice variables We are decently balanced along agricultural practice variables. Our course of action here is similiar to our options with the demographic variables.

Soil Variables We are balanced on the primary soil variables of interest betwen our Tubura farmers and the comparison farmers.

3.4 Baseline balance by district

dist.output <- do.call(rbind, lapply(split(d, d$district), function(x) {
  
  tab <- do.call(rbind, lapply(out.list, function(y) {
    
    out <- t.test(x[,y] ~ x[,"client"], data=x)
    tab <- data.frame(out[[5]][[2]],out[[5]][[1]], out[3])
    tab[,1:2] <- round(tab[,1:2],3)
    names(tab) <- c(names(out[[5]]), "pvalue")
    #tab[,3] <- p.adjust(tab[,3], method="holm")
    #tab[,3] <- ifelse(tab[,3] < 0.001, "< 0.001", round(tab[,3],3))
    #print(tab)
    return(tab)
  }))
  
  return(data.frame(district = unique(x$district), tab))
}))
rownames(dist.output) <- NULL
dist.output$variable <- rep(out.list,length(unique(d$district)))    
# order variables 
dist.output <- dist.output[, c(1, 5, 2:4)]
dist.output$pvalue <- p.adjust(dist.output$pvalue, method="fdr")
dist.output <- dist.output[order(dist.output$pvalue),]
dist.output$pvalue <- ifelse(dist.output$pvalue < 0.001, "< 0.001", round(dist.output$pvalue,3))
colnames(dist.output) <- c("District", "Varible", "Tubura Client","Non-Tubura", "p-value")  

	District	Varible	Tubura Client	Non-Tubura	p-value
162	Karongi	n_season_fert	3.785	0.447	< 0.001
240	Mugonero	n_season_fert	4.313	0.597	< 0.001
435	Rutsiro	n_season_fert	4.200	1.066	< 0.001
318	Nyamasheke	n_season_fert	5.068	1.507	< 0.001
474	Rwamagana	n_season_fert	2.645	1.054	< 0.001
123	Huye	n_season_fert	2.246	0.322	< 0.001
45	Gatsibo_NLWH	n_season_fert	1.517	0.296	< 0.001
461	Rutsiro	own.cows	0.824	0.545	< 0.001
279	Nyamagabe	n_season_fert	4.000	1.554	< 0.001
131	Huye	valley	0.254	0.052	< 0.001
132	Huye	hillside	0.702	0.922	< 0.001
320	Nyamasheke	n_season_lime	0.327	0.041	< 0.001
201	Kayonza	n_season_fert	2.073	0.776	0.001
470	Rwamagana	age	44.245	52.072	0.002
315	Nyamasheke	hhsize	5.701	4.685	0.006
157	Karongi	female	0.466	0.667	0.008
393	Nyaruguru	hhsize	6.217	4.717	0.008
110	Gisagara	own.cows	0.696	0.326	0.008
112	Gisagara	own.chickens	0.370	0.065	0.009
113	Gisagara	own.pigs	0.565	0.217	0.012
119	Huye	age	43.395	49.400	0.017
95	Gisagara	cows	1.043	0.413	0.017
105	Gisagara	m3.Mg	191.739	239.872	0.021
352	Nyanza	female	0.217	0.543	0.022
104	Gisagara	m3.Ca	903.287	1203.933	0.023
357	Nyanza	n_season_fert	1.326	0.152	0.033
41	Gatsibo_NLWH	age	39.897	46.210	0.035
40	Gatsibo_NLWH	female	0.322	0.556	0.036
97	Gisagara	chickens	1.196	0.174	0.036
502	Rwamagana	own.chickens	0.418	0.225	0.036
2	Gatsibo_LWH	age	41.931	49.949	0.046
274	Nyamagabe	female	0.482	0.750	0.052
196	Kayonza	female	0.400	0.672	0.052
149	Huye	own.cows	0.675	0.487	0.054
43	Gatsibo_NLWH	own	0.908	1.000	0.055
242	Mugonero	n_season_lime	0.176	0.023	0.055
469	Rwamagana	female	0.555	0.739	0.055
1	Gatsibo_LWH	female	0.241	0.492	0.06
432	Rutsiro	hhsize	5.624	4.844	0.06
334	Nyamasheke	climbing_bean	0.422	0.267	0.065
305	Nyamagabe	own.cows	0.732	0.482	0.077
337	Nyamasheke	maize	0.082	0.014	0.077
181	Karongi	maize	0.086	0.019	0.078
237	Mugonero	hhsize	5.802	5.054	0.08
310	Nyamagabe	black_soil	0.446	0.214	0.097
62	Gatsibo_NLWH	sorghum	0.276	0.469	0.104
236	Mugonero	age	44.290	49.140	0.107
107	Gisagara	Total.N	0.142	0.155	0.123
134	Huye	cows	1.167	0.774	0.123
391	Nyaruguru	female	0.391	0.652	0.123
324	Nyamasheke	slope	14.129	17.199	0.131
437	Rutsiro	n_season_lime	0.327	0.096	0.131
108	Gisagara	Total.C	1.815	2.005	0.132
319	Nyamasheke	n_season_compost	6.544	5.521	0.14
365	Nyanza	valley	0.283	0.087	0.14
366	Nyanza	hillside	0.717	0.913	0.14
322	Nyamasheke	n_seasons_leg_1	2.327	1.658	0.159
500	Rwamagana	own.cows	0.555	0.396	0.159
256	Mugonero	climbing_bean	0.321	0.194	0.162
396	Nyaruguru	n_season_fert	2.130	1.065	0.162
98	Gisagara	pigs	0.609	0.304	0.17
137	Huye	pigs	0.965	0.748	0.17
81	Gisagara	hhsize	5.348	4.326	0.173
6	Gatsibo_LWH	n_season_fert	2.552	1.475	0.176
128	Huye	n_seasons_leg_2	3.667	4.860	0.176
140	Huye	sorghum	0.237	0.374	0.176
363	Nyanza	slope	7.435	9.239	0.176
426	Nyaruguru	own.sheep	0.109	0.000	0.176
491	Rwamagana	sorghum	0.191	0.324	0.176
152	Huye	own.pigs	0.754	0.617	0.181
450	Rutsiro	sheep	0.467	0.246	0.189
120	Huye	hhsize	5.132	4.574	0.2
276	Nyamagabe	hhsize	5.482	4.696	0.2
333	Nyamasheke	sheep	0.082	0.014	0.211
292	Nyamagabe	chickens	0.857	0.304	0.216
463	Rutsiro	own.chickens	0.291	0.192	0.227
307	Nyamagabe	own.chickens	0.286	0.125	0.231
441	Rutsiro	slope	15.521	13.958	0.235
191	Karongi	own.pigs	0.399	0.289	0.241
207	Kayonza	slope	3.327	1.690	0.241
87	Gisagara	n_season_fallow	0.500	0.087	0.242
89	Gisagara	n_seasons_leg_2	2.370	3.283	0.242
245	Mugonero	n_seasons_leg_2	3.131	4.031	0.246
59	Gatsibo_NLWH	pigs	0.241	0.062	0.253
106	Gisagara	pH	5.759	5.950	0.253
355	Nyanza	own	0.913	1.000	0.256
472	Rwamagana	own	0.909	0.973	0.256
109	Gisagara	ExAl	0.337	0.191	0.273
445	Rutsiro	hilltop	0.097	0.042	0.273
485	Rwamagana	cows	1.082	0.514	0.273
503	Rwamagana	own.pigs	0.282	0.171	0.273
487	Rwamagana	chickens	1.609	0.955	0.28
50	Gatsibo_NLWH	n_seasons_leg_2	3.046	3.889	0.281
84	Gisagara	n_season_fert	1.087	0.348	0.281
348	Nyamasheke	own.sheep	0.054	0.014	0.288
362	Nyanza	n_seasons_leg_2	2.130	3.478	0.289
164	Karongi	n_season_lime	0.098	0.013	0.29
203	Kayonza	n_season_lime	0.527	0.707	0.293
370	Nyanza	chickens	2.348	1.174	0.293
477	Rwamagana	n_season_fallow	1.218	0.712	0.293
411	Nyaruguru	sheep	0.196	0.000	0.295
96	Gisagara	goats	1.522	0.935	0.296
275	Nyamagabe	age	43.768	48.875	0.296
172	Karongi	hilltop	0.067	0.025	0.338
92	Gisagara	valley	0.283	0.130	0.338
294	Nyamagabe	sheep	0.286	0.089	0.338
407	Nyaruguru	cows	1.196	0.891	0.338
419	Nyaruguru	Total.N	0.161	0.170	0.338
347	Nyamasheke	own.pigs	0.408	0.308	0.342
146	Huye	Total.N	0.142	0.148	0.352
173	Karongi	cows	1.387	1.119	0.374
55	Gatsibo_NLWH	hilltop	0.011	0.062	0.381
142	Huye	maize	0.061	0.017	0.381
216	Kayonza	sheep	0.036	0.172	0.381
233	Kayonza	red_soil	0.200	0.086	0.381
382	Nyanza	ExAl	0.117	0.175	0.381
353	Nyanza	age	44.717	50.152	0.397
422	Nyaruguru	own.cows	0.826	0.674	0.397
79	Gisagara	female	0.435	0.609	0.403
425	Nyaruguru	own.pigs	0.609	0.435	0.403
158	Karongi	age	44.681	47.384	0.404
225	Kayonza	Total.C	2.395	2.552	0.408
410	Nyaruguru	pigs	0.783	0.457	0.408
495	Rwamagana	m3.Mg	270.948	289.488	0.409
266	Mugonero	own.cows	0.786	0.698	0.413
197	Kayonza	age	42.345	46.500	0.415
312	Nyamagabe	sandy_soil	0.089	0.196	0.421
436	Rutsiro	n_season_compost	8.176	7.593	0.441
46	Gatsibo_NLWH	n_season_compost	2.632	3.420	0.442
33	Gatsibo_LWH	own.goats	0.552	0.407	0.455
221	Kayonza	m3.Ca	1677.807	1878.700	0.458
375	Nyanza	bush_bean	0.261	0.413	0.473
180	Karongi	bush_bean	0.080	0.132	0.481
446	Rutsiro	cows	1.442	1.048	0.485
424	Nyaruguru	own.chickens	0.283	0.152	0.488
260	Mugonero	m3.Ca	551.499	603.770	0.49
390	Nyanza	sandy_soil	0.304	0.457	0.493
438	Rutsiro	n_season_fallow	0.473	0.287	0.511
323	Nyamasheke	n_seasons_leg_2	2.137	2.681	0.515
32	Gatsibo_LWH	own.cows	0.500	0.373	0.526
74	Gatsibo_NLWH	own.pigs	0.126	0.062	0.526
102	Gisagara	bush_bean	0.217	0.348	0.526
103	Gisagara	maize	0.043	0.000	0.526
118	Huye	female	0.500	0.591	0.526
153	Huye	own.sheep	0.053	0.017	0.526
199	Kayonza	own	0.964	0.897	0.526
200	Kayonza	field.size	664.236	423.776	0.526
231	Kayonza	own.sheep	0.036	0.103	0.526
250	Mugonero	hilltop	0.015	0.000	0.526
283	Nyamagabe	n_seasons_leg_1	1.625	1.125	0.526
293	Nyamagabe	pigs	1.000	0.804	0.526
297	Nyamagabe	bush_bean	0.000	0.036	0.526
309	Nyamagabe	own.sheep	0.179	0.089	0.526
379	Nyanza	pH	6.104	5.995	0.526
398	Nyaruguru	n_season_lime	0.043	0.000	0.526
418	Nyaruguru	pH	5.320	5.205	0.526
421	Nyaruguru	ExAl	0.717	0.882	0.526
430	Rutsiro	female	0.545	0.623	0.526
488	Rwamagana	pigs	0.518	0.315	0.526
507	Rwamagana	sandy_soil	0.018	0.000	0.526
47	Gatsibo_NLWH	n_season_lime	0.069	0.025	0.533
349	Nyamasheke	black_soil	0.531	0.610	0.533
52	Gatsibo_NLWH	alt	1017.040	1129.172	0.534
58	Gatsibo_NLWH	chickens	1.678	1.000	0.534
117	Gisagara	sandy_soil	0.239	0.370	0.534
174	Karongi	goats	0.890	1.069	0.534
224	Kayonza	Total.N	0.180	0.187	0.534
498	Rwamagana	Total.C	2.174	2.236	0.534
17	Gatsibo_LWH	cows	0.759	0.542	0.534
371	Nyanza	pigs	0.087	0.217	0.534
386	Nyanza	own.pigs	0.065	0.152	0.534
394	Nyaruguru	own	0.848	0.935	0.534
465	Rutsiro	own.sheep	0.212	0.156	0.535
220	Kayonza	maize	0.018	0.069	0.535
311	Nyamagabe	red_soil	0.464	0.589	0.536
457	Rutsiro	pH	5.269	5.339	0.536
93	Gisagara	hillside	0.587	0.717	0.543
163	Karongi	n_season_compost	7.074	6.553	0.543
76	Gatsibo_NLWH	black_soil	0.609	0.704	0.554
214	Kayonza	chickens	0.982	0.466	0.561
455	Rutsiro	m3.Ca	628.725	682.389	0.566
314	Nyamasheke	age	45.898	48.123	0.568
383	Nyanza	own.cows	0.652	0.522	0.568
111	Gisagara	own.goats	0.630	0.500	0.571
160	Karongi	own	0.957	0.981	0.571
222	Kayonza	m3.Mg	331.923	356.976	0.571
262	Mugonero	pH	5.180	5.243	0.571
265	Mugonero	ExAl	0.925	0.823	0.571
344	Nyamasheke	own.cows	0.639	0.568	0.571
377	Nyanza	m3.Ca	1086.559	992.934	0.571
11	Gatsibo_LWH	n_seasons_leg_2	1.983	2.690	0.574
18	Gatsibo_LWH	goats	1.345	0.983	0.581
127	Huye	n_seasons_leg_1	1.386	1.043	0.581
209	Kayonza	valley	0.545	0.431	0.581
378	Nyanza	m3.Mg	233.315	213.344	0.581
19	Gatsibo_LWH	chickens	1.052	0.678	0.583
136	Huye	chickens	1.772	1.348	0.583
442	Rutsiro	alt	1945.122	1894.425	0.583
479	Rwamagana	n_seasons_leg_2	1.349	1.685	0.583
78	Gatsibo_NLWH	sandy_soil	0.287	0.210	0.587
159	Karongi	hhsize	5.172	4.912	0.587
171	Karongi	hillside	0.785	0.836	0.587
218	Kayonza	sorghum	0.182	0.103	0.587
298	Nyamagabe	maize	0.089	0.036	0.587
389	Nyanza	red_soil	0.109	0.043	0.587
417	Nyaruguru	m3.Mg	152.908	139.740	0.587
420	Nyaruguru	Total.C	2.145	2.237	0.587
447	Rutsiro	goats	0.515	0.659	0.587
459	Rutsiro	Total.C	2.247	2.309	0.587
480	Rwamagana	slope	3.682	4.234	0.587
462	Rutsiro	own.goats	0.273	0.329	0.619
340	Nyamasheke	pH	5.176	5.120	0.623
219	Kayonza	bush_bean	0.218	0.138	0.632
234	Kayonza	sandy_soil	0.036	0.086	0.632
64	Gatsibo_NLWH	maize	0.092	0.049	0.648
268	Mugonero	own.chickens	0.382	0.318	0.648
460	Rutsiro	ExAl	0.814	0.736	0.648
37	Gatsibo_LWH	black_soil	0.172	0.254	0.648
489	Rwamagana	sheep	0.509	0.360	0.649
3	Gatsibo_LWH	hhsize	5.362	5.763	0.651
4	Gatsibo_LWH	own	0.983	1.000	0.651
38	Gatsibo_LWH	red_soil	0.638	0.542	0.651
57	Gatsibo_NLWH	goats	1.218	1.654	0.651
71	Gatsibo_NLWH	own.cows	0.563	0.481	0.651
88	Gisagara	n_seasons_leg_1	2.304	2.870	0.651
115	Gisagara	black_soil	0.587	0.478	0.651
121	Huye	own	0.974	0.991	0.651
125	Huye	n_season_lime	0.035	0.113	0.651
135	Huye	goats	0.965	0.809	0.651
154	Huye	black_soil	0.518	0.583	0.651
186	Karongi	Total.C	1.887	1.839	0.651
213	Kayonza	goats	1.727	1.397	0.651
217	Kayonza	climbing_bean	0.000	0.017	0.651
235	Mugonero	female	0.656	0.713	0.651
238	Mugonero	own	0.947	0.915	0.651
259	Mugonero	maize	0.000	0.008	0.651
290	Nyamagabe	cows	0.911	0.714	0.651
343	Nyamasheke	ExAl	0.983	1.062	0.651
350	Nyamasheke	red_soil	0.320	0.267	0.651
368	Nyanza	cows	1.174	0.891	0.651
376	Nyanza	maize	0.130	0.065	0.651
385	Nyanza	own.chickens	0.543	0.435	0.651
395	Nyaruguru	field.size	570.022	495.239	0.651
414	Nyaruguru	bush_bean	0.065	0.022	0.651
415	Nyaruguru	maize	0.022	0.000	0.651
434	Rutsiro	field.size	584.327	430.156	0.651
449	Rutsiro	pigs	0.315	0.240	0.651
493	Rwamagana	maize	0.027	0.009	0.651
494	Rwamagana	m3.Ca	1303.126	1380.293	0.651
497	Rwamagana	Total.N	0.175	0.179	0.651
53	Gatsibo_NLWH	valley	0.759	0.691	0.654
178	Karongi	climbing_bean	0.405	0.352	0.654
251	Mugonero	cows	1.328	1.171	0.654
374	Nyanza	sorghum	0.283	0.196	0.654
443	Rutsiro	valley	0.152	0.192	0.654
211	Kayonza	hilltop	0.109	0.172	0.655
351	Nyamasheke	sandy_soil	0.116	0.082	0.657
56	Gatsibo_NLWH	cows	0.885	0.741	0.658
325	Nyamasheke	alt	1678.150	1704.388	0.658
456	Rutsiro	m3.Mg	165.214	175.386	0.658
228	Kayonza	own.goats	0.673	0.586	0.659
8	Gatsibo_LWH	n_season_lime	0.500	0.322	0.66
83	Gisagara	field.size	493.957	559.261	0.669
141	Huye	bush_bean	0.193	0.243	0.669
302	Nyamagabe	Total.N	0.169	0.165	0.669
409	Nyaruguru	chickens	0.652	0.391	0.669
73	Gatsibo_NLWH	own.chickens	0.310	0.247	0.673
454	Rutsiro	maize	0.176	0.216	0.673
5	Gatsibo_LWH	field.size	776.345	916.441	0.697
91	Gisagara	alt	1036.055	923.525	0.697
303	Nyamagabe	Total.C	2.345	2.274	0.697
364	Nyanza	alt	1505.653	1529.554	0.697
239	Mugonero	field.size	441.359	393.124	0.698
223	Kayonza	pH	6.132	6.204	0.699
458	Rutsiro	Total.N	0.157	0.159	0.699
504	Rwamagana	own.sheep	0.227	0.180	0.699
194	Karongi	red_soil	0.325	0.371	0.699
34	Gatsibo_LWH	own.chickens	0.345	0.271	0.702
150	Huye	own.goats	0.518	0.461	0.702
48	Gatsibo_NLWH	n_season_fallow	0.736	0.543	0.703
332	Nyamasheke	pigs	0.776	0.582	0.703
423	Nyaruguru	own.goats	0.348	0.435	0.704
384	Nyanza	own.goats	0.630	0.543	0.706
405	Nyaruguru	hillside	0.957	0.913	0.706
406	Nyaruguru	hilltop	0.043	0.087	0.706
101	Gisagara	sorghum	0.435	0.522	0.709
291	Nyamagabe	goats	0.750	0.607	0.709
313	Nyamasheke	female	0.646	0.692	0.709
167	Karongi	n_seasons_leg_2	3.216	3.548	0.71
308	Nyamagabe	own.pigs	0.732	0.661	0.714
63	Gatsibo_NLWH	bush_bean	0.172	0.222	0.715
126	Huye	n_season_fallow	0.675	0.487	0.715
317	Nyamasheke	field.size	787.527	580.130	0.715
354	Nyanza	hhsize	5.261	4.891	0.715
416	Nyaruguru	m3.Ca	660.083	610.274	0.715
254	Mugonero	pigs	0.412	0.333	0.719
227	Kayonza	own.cows	0.382	0.310	0.72
189	Karongi	own.goats	0.497	0.541	0.721
269	Mugonero	own.pigs	0.244	0.287	0.73
403	Nyaruguru	alt	1791.131	1814.738	0.73
272	Mugonero	red_soil	0.252	0.295	0.73
278	Nyamagabe	field.size	310.929	282.643	0.73
80	Gisagara	age	44.761	46.848	0.736
270	Mugonero	own.sheep	0.122	0.093	0.736
492	Rwamagana	bush_bean	0.364	0.315	0.736
155	Huye	red_soil	0.175	0.139	0.737
261	Mugonero	m3.Mg	160.377	170.165	0.74
10	Gatsibo_LWH	n_seasons_leg_1	4.345	4.678	0.747
23	Gatsibo_LWH	sorghum	0.276	0.339	0.747
26	Gatsibo_LWH	m3.Ca	1226.445	1149.970	0.747
133	Huye	hilltop	0.044	0.026	0.747
161	Karongi	field.size	490.788	446.088	0.747
176	Karongi	pigs	0.528	0.447	0.747
429	Nyaruguru	sandy_soil	0.109	0.065	0.747
187	Karongi	ExAl	0.647	0.597	0.752
215	Kayonza	pigs	0.309	0.224	0.753
138	Huye	sheep	0.061	0.035	0.754
496	Rwamagana	pH	5.817	5.862	0.764
44	Gatsibo_NLWH	field.size	1306.069	1459.901	0.772
124	Huye	n_season_compost	5.211	4.878	0.772
24	Gatsibo_LWH	bush_bean	0.397	0.339	0.787
60	Gatsibo_NLWH	sheep	0.046	0.025	0.787
72	Gatsibo_NLWH	own.goats	0.517	0.568	0.787
183	Karongi	m3.Mg	254.160	267.784	0.787
185	Karongi	Total.N	0.144	0.142	0.787
188	Karongi	own.cows	0.767	0.736	0.787
205	Kayonza	n_seasons_leg_1	1.982	2.155	0.787
341	Nyamasheke	Total.N	0.165	0.167	0.787
360	Nyanza	n_season_fallow	1.130	0.870	0.787
427	Nyaruguru	black_soil	0.630	0.696	0.787
433	Rutsiro	own	0.915	0.934	0.787
20	Gatsibo_LWH	pigs	0.138	0.203	0.79
85	Gisagara	n_season_compost	4.326	3.957	0.79
229	Kayonza	own.chickens	0.182	0.138	0.79
184	Karongi	pH	5.442	5.477	0.79
388	Nyanza	black_soil	0.565	0.500	0.79
467	Rutsiro	red_soil	0.315	0.347	0.79
475	Rwamagana	n_season_compost	4.018	3.748	0.79
271	Mugonero	black_soil	0.527	0.488	0.791
230	Kayonza	own.pigs	0.218	0.172	0.794
244	Mugonero	n_seasons_leg_1	1.336	1.186	0.794
408	Nyaruguru	goats	0.609	0.739	0.794
482	Rwamagana	valley	0.409	0.369	0.794
483	Rwamagana	hillside	0.527	0.568	0.795
36	Gatsibo_LWH	own.sheep	0.034	0.017	0.796
42	Gatsibo_NLWH	hhsize	5.333	5.160	0.796
206	Kayonza	n_seasons_leg_2	1.091	1.259	0.796
277	Nyamagabe	own	0.964	0.982	0.796
372	Nyanza	sheep	0.022	0.043	0.796
373	Nyanza	climbing_bean	0.043	0.022	0.796
387	Nyanza	own.sheep	0.022	0.043	0.796
399	Nyaruguru	n_season_fallow	0.370	0.261	0.796
486	Rwamagana	goats	1.745	1.568	0.796
82	Gisagara	own	0.826	0.870	0.797
9	Gatsibo_LWH	n_season_fallow	0.397	0.288	0.801
280	Nyamagabe	n_season_compost	7.857	7.518	0.801
306	Nyamagabe	own.goats	0.464	0.411	0.801
431	Rutsiro	age	43.685	44.617	0.801
210	Kayonza	hillside	0.345	0.397	0.802
90	Gisagara	slope	8.391	7.739	0.807
339	Nyamasheke	m3.Mg	156.068	149.714	0.807
281	Nyamagabe	n_season_lime	0.500	0.393	0.808
304	Nyamagabe	ExAl	1.111	1.170	0.814
143	Huye	m3.Ca	850.319	872.728	0.821
144	Huye	m3.Mg	187.671	191.878	0.821
286	Nyamagabe	alt	2124.389	2109.814	0.824
66	Gatsibo_NLWH	m3.Mg	246.996	255.132	0.824
267	Mugonero	own.goats	0.542	0.574	0.824
122	Huye	field.size	599.096	567.922	0.828
284	Nyamagabe	n_seasons_leg_2	3.161	2.929	0.828
299	Nyamagabe	m3.Ca	458.692	440.797	0.828
412	Nyaruguru	climbing_bean	0.196	0.239	0.828
468	Rutsiro	sandy_soil	0.176	0.156	0.833
77	Gatsibo_NLWH	red_soil	0.034	0.049	0.842
338	Nyamasheke	m3.Ca	614.196	590.190	0.842
505	Rwamagana	black_soil	0.655	0.685	0.842
192	Karongi	own.sheep	0.080	0.094	0.842
246	Mugonero	slope	15.939	16.473	0.842
316	Nyamasheke	own	0.932	0.945	0.842
397	Nyaruguru	n_season_compost	5.696	6.043	0.842
501	Rwamagana	own.goats	0.536	0.568	0.842
151	Huye	own.chickens	0.421	0.391	0.845
296	Nyamagabe	sorghum	0.036	0.054	0.845
356	Nyanza	field.size	853.359	766.217	0.845
253	Mugonero	chickens	1.076	0.961	0.847
400	Nyaruguru	n_seasons_leg_1	2.043	2.348	0.847
361	Nyanza	n_seasons_leg_1	2.543	2.283	0.858
466	Rutsiro	black_soil	0.461	0.437	0.858
39	Gatsibo_LWH	sandy_soil	0.172	0.203	0.859
190	Karongi	own.chickens	0.411	0.434	0.862
473	Rwamagana	field.size	892.855	842.054	0.862
478	Rwamagana	n_seasons_leg_1	1.682	1.568	0.862
336	Nyamasheke	bush_bean	0.088	0.075	0.867
175	Karongi	chickens	1.503	1.377	0.87
476	Rwamagana	n_season_lime	0.355	0.324	0.87
490	Rwamagana	climbing_bean	0.036	0.027	0.872
75	Gatsibo_NLWH	own.sheep	0.034	0.025	0.89
330	Nyamasheke	goats	0.918	0.863	0.891
177	Karongi	sheep	0.166	0.195	0.903
195	Karongi	sandy_soil	0.025	0.019	0.906
13	Gatsibo_LWH	alt	1094.532	1126.650	0.911
145	Huye	pH	5.681	5.664	0.911
241	Mugonero	n_season_compost	6.229	6.070	0.911
464	Rutsiro	own.pigs	0.194	0.180	0.911
321	Nyamasheke	n_season_fallow	0.599	0.664	0.911
35	Gatsibo_LWH	own.pigs	0.121	0.102	0.916
444	Rutsiro	hillside	0.752	0.766	0.919
28	Gatsibo_LWH	pH	6.139	6.115	0.925
29	Gatsibo_LWH	Total.N	0.148	0.147	0.925
30	Gatsibo_LWH	Total.C	1.882	1.856	0.925
31	Gatsibo_LWH	ExAl	0.151	0.164	0.925
94	Gisagara	hilltop	0.130	0.152	0.925
116	Gisagara	red_soil	0.152	0.130	0.925
129	Huye	slope	6.658	6.817	0.925
193	Karongi	black_soil	0.583	0.566	0.925
249	Mugonero	hillside	0.855	0.868	0.925
287	Nyamagabe	valley	0.125	0.107	0.925
130	Huye	alt	1647.892	1655.403	0.928
342	Nyamasheke	Total.C	2.307	2.289	0.928
401	Nyaruguru	n_seasons_leg_2	3.178	3.370	0.928
481	Rwamagana	alt	982.728	961.878	0.928
54	Gatsibo_NLWH	hillside	0.230	0.247	0.931
67	Gatsibo_NLWH	pH	6.037	6.053	0.931
156	Huye	sandy_soil	0.263	0.278	0.931
182	Karongi	m3.Ca	781.863	796.295	0.931
295	Nyamagabe	climbing_bean	0.161	0.143	0.931
301	Nyamagabe	pH	4.999	5.013	0.931
288	Nyamagabe	hillside	0.821	0.839	0.935
439	Rutsiro	n_seasons_leg_1	3.388	3.479	0.935
170	Karongi	valley	0.147	0.138	0.939
212	Kayonza	cows	0.727	0.828	0.939
243	Mugonero	n_season_fallow	0.863	0.915	0.939
257	Mugonero	sorghum	0.061	0.054	0.939
380	Nyanza	Total.N	0.134	0.135	0.939
428	Nyaruguru	red_soil	0.261	0.239	0.939
27	Gatsibo_LWH	m3.Mg	238.242	234.655	0.942
139	Huye	climbing_bean	0.088	0.096	0.942
226	Kayonza	ExAl	0.111	0.117	0.942
285	Nyamagabe	slope	12.571	12.839	0.942
300	Nyamagabe	m3.Mg	109.343	111.050	0.942
331	Nyamasheke	chickens	0.599	0.562	0.942
499	Rwamagana	ExAl	0.257	0.249	0.942
70	Gatsibo_NLWH	ExAl	0.182	0.175	0.945
282	Nyamagabe	n_season_fallow	0.268	0.304	0.945
381	Nyanza	Total.C	1.766	1.782	0.945
148	Huye	ExAl	0.307	0.315	0.948
232	Kayonza	black_soil	0.691	0.707	0.948
471	Rwamagana	hhsize	4.927	4.982	0.948
506	Rwamagana	red_soil	0.327	0.315	0.948
448	Rutsiro	chickens	0.879	0.958	0.949
147	Huye	Total.C	1.917	1.925	0.95
65	Gatsibo_NLWH	m3.Ca	1236.151	1251.746	0.951
198	Kayonza	hhsize	5.091	5.155	0.951
21	Gatsibo_LWH	sheep	0.103	0.085	0.955
451	Rutsiro	climbing_bean	0.327	0.335	0.959
7	Gatsibo_LWH	n_season_compost	5.672	5.780	0.965
247	Mugonero	alt	1796.175	1788.315	0.965
326	Nyamasheke	valley	0.150	0.144	0.966
51	Gatsibo_NLWH	slope	4.011	4.099	0.966
69	Gatsibo_NLWH	Total.C	2.002	1.990	0.966
204	Kayonza	n_season_fallow	0.455	0.483	0.967
440	Rutsiro	n_seasons_leg_2	2.036	2.078	0.967
327	Nyamasheke	hillside	0.823	0.829	0.967
14	Gatsibo_LWH	valley	0.603	0.593	0.976
15	Gatsibo_LWH	hillside	0.397	0.407	0.976
255	Mugonero	sheep	0.221	0.233	0.98
68	Gatsibo_NLWH	Total.N	0.157	0.158	0.981
179	Karongi	sorghum	0.147	0.151	0.983
252	Mugonero	goats	1.321	1.302	0.983
358	Nyanza	n_season_compost	3.261	3.326	0.984
369	Nyanza	goats	1.261	1.283	0.994
263	Mugonero	Total.N	0.153	0.153	0.997
202	Kayonza	n_season_compost	3.564	3.534	0.999
208	Kayonza	alt	1196.807	1202.646	0.999
258	Mugonero	bush_bean	0.176	0.178	0.999
273	Mugonero	sandy_soil	0.176	0.178	0.999
248	Mugonero	valley	0.130	0.132	1
345	Nyamasheke	own.goats	0.442	0.445	1
392	Nyaruguru	age	48.457	48.304	1
12	Gatsibo_LWH	slope	4.259	4.220	1
49	Gatsibo_NLWH	n_seasons_leg_1	1.540	1.556	1
86	Gisagara	n_season_lime	0.022	0.022	1
99	Gisagara	sheep	0.022	0.022	1
100	Gisagara	climbing_bean	0.022	0.022	1
114	Gisagara	own.sheep	0.022	0.022	1
165	Karongi	n_season_fallow	0.834	0.843	1
166	Karongi	n_seasons_leg_1	2.485	2.497	1
168	Karongi	slope	11.975	11.937	1
169	Karongi	alt	1806.841	1806.879	1
264	Mugonero	Total.C	2.217	2.218	1
289	Nyamagabe	hilltop	0.054	0.054	1
328	Nyamasheke	hilltop	0.027	0.027	1
329	Nyamasheke	cows	1.245	1.240	1
346	Nyamasheke	own.chickens	0.218	0.219	1
402	Nyaruguru	slope	9.478	9.457	1
413	Nyaruguru	sorghum	0.304	0.304	1
453	Rutsiro	bush_bean	0.012	0.012	1
484	Rwamagana	hilltop	0.064	0.063	1
16	Gatsibo_LWH	hilltop	0.000	0.000	NA
22	Gatsibo_LWH	climbing_bean	0.000	0.000	NA
25	Gatsibo_LWH	maize	0.000	0.000	NA
61	Gatsibo_NLWH	climbing_bean	0.000	0.000	NA
335	Nyamasheke	sorghum	0.000	0.000	NA
359	Nyanza	n_season_lime	0.000	0.000	NA
367	Nyanza	hilltop	0.000	0.000	NA
404	Nyaruguru	valley	0.000	0.000	NA
452	Rutsiro	sorghum	0.000	0.000	NA

3.5 District balance interpretation

Demographic variables interpretation here.

Agricultural practice variables interpretation here

Soil Variables interpretation here

3.6 Baseline balance by Tubura tenure

Look at farmers by duration of tenure farming with Tubura. We want to understand, at least with an initial naive baseline sense, what is the cumulative effect of Tubura practices on soil health outcomes?

We will look only at current Tubura farmers and compare first year farmers to farmers with more experience with Tubura.

oafOnly <- d[which(d$client==1 & d$total.seasons>=1),]
nTenure <- oafOnly %>% group_by(total.seasons) %>% 
  dplyr::summarize(
    n = n()
  ) %>% ungroup() %>% as.data.frame()
nTenure$val <- paste(nTenure$total.seasons, " (", "n = ", nTenure$n, ")", sep = "")
for(i in 1:length(soilVars)){
print(
  ggplot(oafOnly, aes(x=as.factor(total.seasons), y=oafOnly[,soilVars[i]])) + 
    geom_boxplot() +
    scale_x_discrete(labels=nTenure$val) + 
    theme(legend.position = "bottom", axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(x="Tubura Tenure", y=soilVars[i], title = paste("RW baseline soil by tenure - ", soilVars[i], sep = ""))
  )  
}

3.7 Tenure summaries

tenureSum <- aggregate(oafOnly[,out.list], by=list(oafOnly$total.seasons), function(x){
  round(mean(x, na.rm=T),2)
})
tenureSum <- as.data.frame(t(tenureSum))
colnames(tenureSum) <- c(paste(seq(1,14,1), " seas.", sep = ""))

print(kable(tenureSum))

	1 seas.	2 seas.	3 seas.	4 seas.	5 seas.	6 seas.	7 seas.	8 seas.	9 seas.	10 seas.	11 seas.	12 seas.	13 seas.	14 seas.
Group.1	1.00	2.00	3.00	4.00	5.00	6.00	7.00	8.00	9.00	10.00	11.00	12.00	13.00	14.00
female	0.42	0.48	0.60	0.49	0.59	0.59	0.52	0.50	0.65	0.62	0.75	0.53	0.55	0.48
age	42.63	42.68	44.33	45.58	43.28	43.96	44.19	50.91	51.85	42.76	48.25	48.97	52.73	49.52
hhsize	5.24	5.17	5.08	5.21	6.07	5.58	6.19	5.66	6.65	6.07	6.58	5.53	4.55	6.12
own	0.93	0.95	0.92	0.95	0.98	0.97	1.00	0.93	0.85	1.00	0.92	0.89	0.91	0.95
field.size	606.74	630.98	630.65	551.74	482.02	484.79	412.52	592.95	330.00	614.79	617.83	440.76	280.82	679.71
n_season_fert	1.76	2.32	2.63	3.28	4.67	4.90	6.07	5.59	6.00	5.83	7.25	7.18	7.18	6.79
n_season_compost	5.35	5.37	4.81	6.37	6.87	7.13	8.30	7.43	7.30	7.52	8.67	8.24	8.27	7.81
n_season_lime	0.16	0.27	0.36	0.31	0.22	0.26	0.07	0.07	0.20	0.28	0.08	0.16	0.09	0.48
n_season_fallow	0.65	0.50	0.77	0.90	0.65	0.66	0.67	0.77	1.10	0.21	0.42	0.21	1.82	0.67
n_seasons_leg_1	1.97	2.07	2.10	2.44	2.96	2.89	1.96	2.91	1.60	1.79	2.92	3.53	2.09	3.02
n_seasons_leg_2	2.44	2.54	2.03	2.32	2.67	2.62	4.93	3.50	2.75	2.59	3.92	1.49	2.82	2.12
slope	9.05	8.66	9.65	8.80	9.57	11.84	14.07	11.77	13.60	12.31	14.17	11.26	14.91	14.38
alt	1531.01	1491.56	1550.41	1518.25	1654.93	1582.27	1758.60	1639.74	1812.98	1764.71	1788.95	1640.20	1702.35	1633.76
valley	0.31	0.32	0.27	0.25	0.30	0.22	0.26	0.23	0.00	0.07	0.33	0.18	0.09	0.21
hillside	0.62	0.61	0.71	0.72	0.67	0.74	0.70	0.66	0.95	0.90	0.67	0.79	0.82	0.76
hilltop	0.07	0.06	0.03	0.03	0.02	0.05	0.04	0.11	0.05	0.03	0.00	0.03	0.09	0.02
cows	0.99	1.00	0.94	1.49	1.91	1.28	1.30	1.09	3.25	1.31	1.75	0.97	1.00	1.07
goats	0.97	1.12	1.00	1.28	1.39	1.00	0.78	1.16	1.20	0.90	0.83	1.11	0.45	1.26
chickens	1.07	1.22	1.08	1.20	1.28	0.87	1.78	1.61	1.20	1.45	3.58	1.76	0.45	1.00
pigs	0.55	0.52	0.58	0.51	0.37	0.47	0.48	0.55	0.45	0.55	0.25	0.84	0.64	0.83
sheep	0.20	0.21	0.29	0.34	0.17	0.16	0.44	0.23	0.15	0.41	0.00	0.05	0.18	0.17
climbing_bean	0.14	0.17	0.21	0.17	0.20	0.32	0.52	0.36	0.35	0.45	0.33	0.39	0.55	0.38
sorghum	0.18	0.17	0.12	0.14	0.11	0.06	0.04	0.02	0.05	0.10	0.17	0.05	0.09	0.00
bush_bean	0.14	0.17	0.17	0.17	0.15	0.19	0.04	0.18	0.15	0.03	0.08	0.11	0.18	0.14
maize	0.06	0.07	0.06	0.08	0.09	0.08	0.11	0.11	0.10	0.10	0.08	0.03	0.00	0.05
m3.Ca	890.32	911.95	871.92	918.19	913.35	840.76	771.94	755.95	592.43	582.74	663.89	848.28	592.16	618.11
m3.Mg	195.33	201.53	213.24	210.61	218.70	214.24	210.69	222.18	156.45	172.28	198.05	275.01	151.55	166.45
pH	5.57	5.57	5.57	5.53	5.64	5.52	5.43	5.46	5.19	5.22	5.37	5.57	5.29	5.25
Total.N	0.16	0.16	0.16	0.16	0.15	0.15	0.15	0.14	0.16	0.16	0.14	0.15	0.15	0.16
Total.C	2.11	2.15	2.11	2.15	2.13	2.07	2.05	1.93	2.38	2.16	1.88	1.97	2.00	2.22
ExAl	0.59	0.57	0.62	0.54	0.46	0.61	0.60	0.56	0.91	0.75	0.70	0.48	0.64	0.90
own.cows	0.69	0.63	0.58	0.75	0.72	0.69	0.81	0.70	0.90	0.90	0.75	0.63	0.64	0.81
own.goats	0.45	0.52	0.47	0.46	0.54	0.50	0.48	0.55	0.50	0.48	0.50	0.53	0.27	0.52
own.chickens	0.30	0.35	0.33	0.30	0.33	0.30	0.56	0.43	0.45	0.45	0.58	0.42	0.45	0.33
own.pigs	0.42	0.36	0.32	0.31	0.26	0.33	0.41	0.36	0.35	0.34	0.17	0.39	0.36	0.40
own.sheep	0.13	0.09	0.14	0.17	0.11	0.08	0.26	0.14	0.10	0.14	0.00	0.03	0.09	0.07
black_soil	0.54	0.54	0.55	0.57	0.50	0.43	0.52	0.45	0.55	0.48	0.25	0.61	0.45	0.62
red_soil	0.29	0.25	0.29	0.25	0.28	0.37	0.26	0.30	0.30	0.34	0.42	0.24	0.36	0.26
sandy_soil	0.13	0.19	0.14	0.15	0.17	0.18	0.07	0.18	0.05	0.07	0.17	0.13	0.18	0.10

3.8 Tenure balance

We’re defining Tubura tenure as having 3 or more seasons of experience farming with Tubura. We draw the line at 3 seasons as three seasons of fertilizer use is approximately when we’d expect fertilizer to start to have a detrimental effect on soil health.

oafOnly$tenured <- ifelse(oafOnly$total.seasons>=3,1,0)
tenure <- do.call(rbind, lapply(out.list, function(x) {
  
  out <- t.test(oafOnly[,x] ~ oafOnly[,"tenured"], data=oafOnly)
  tab <- data.frame(out[[5]][[1]], out[[5]][[2]], out[3])
  tab[,1:2] <- round(tab[,1:2],3)
  names(tab) <- c(names(out[[5]]), "pvalue")
  return(tab)
}))
# use p.adjust with bonferroni correction
tenure$pvalue <- p.adjust(tenure$pvalue, method="fdr")
rownames(tenure) <- out.list
tenure <- tenure[order(tenure$pvalue),]
tenure$pvalue <- ifelse(tenure[, 3] < 0.001, "< 0.001", round(tenure[, 3], 3)) 
colnames(tenure) <- c("Non-Tubura", "Tubura Client", "p-value") 

print(kable(tenure))

	Non-Tubura	Tubura Client	p-value
n_season_fert	2.097	4.774	< 0.001
n_season_compost	5.359	6.899	< 0.001
climbing_bean	0.159	0.295	< 0.001
age	42.657	46.057	< 0.001
slope	8.814	11.103	< 0.001
hillside	0.616	0.736	< 0.001
sorghum	0.173	0.086	< 0.001
n_seasons_leg_1	2.030	2.583	0.006
valley	0.319	0.225	0.007
hhsize	5.200	5.597	0.011
cows	0.992	1.358	0.011
female	0.457	0.552	0.013
alt	1507.236	1611.057	0.017
m3.Ca	903.359	813.631	0.037
pH	5.570	5.478	0.039
Total.N	0.159	0.154	0.047
own.cows	0.657	0.719	0.105
field.size	621.349	535.512	0.127
hilltop	0.065	0.039	0.193
own.pigs	0.386	0.334	0.204
n_season_fallow	0.559	0.719	0.227
m3.Mg	199.065	208.991	0.362
Total.C	2.136	2.102	0.505
sandy_soil	0.168	0.142	0.505
ExAl	0.574	0.613	0.52
own.sheep	0.103	0.122	0.52
red_soil	0.268	0.295	0.52
own.chickens	0.332	0.360	0.529
sheep	0.203	0.242	0.603
maize	0.065	0.077	0.607
n_season_lime	0.227	0.259	0.61
black_soil	0.541	0.518	0.61
chickens	1.162	1.261	0.624
goats	1.062	1.111	0.717
own	0.941	0.947	0.762
n_seasons_leg_2	2.500	2.575	0.764
bush_bean	0.159	0.151	0.764
pigs	0.532	0.545	0.869
own.goats	0.492	0.492	0.991

3.9 Tenured v. new balance interpretation

Demographic variables We are well balanced along demographic variables.

Agricultural practice variables Not surprisingly, Tubura farmers have more cumulative years of fertilizer use than current non-Tubura farmers. While that difference is signficant, it is realistically only a single season of fertilizer use different.

Interestingly, non-Tubura farmers reported using more lime than current Tubura farmers. This

Soil Variables Soil pH, calcium and magnesium levels are lower for tenured Tubura farmers. This is consistent with the hypothesis that increaesd fertilizer use leads to an increaese in soil acidity.

3.10 Analysis of agronomic practices contributing to soil health outcomes

Here’s where we’ll look at the contribution of fertilizer, lime and cultivation practices on soil health outcomes. This analysis will be come richer as we gain longitudinal measures. I caution that we cannot treat these relationships as causal. The direction of causality is not clearly delineated in the data or the study design. However, we can identify meaningful connections between practices and outcomes through this analysis to generate new hypotheses for field testing.

I’m going to start with behaviors by sections and then move to a more comprehensive model including multiple practices. All models will include controls for site to account for local variation and field officer behavior.

4 Analysis

4.1 Agronomic Behaviors

Check for multicollinearity before adding number of seasons of agronomic inputs on the same side of the regression.

suppressMessages(library(stargazer))
inputUse <- c("n_season_fert","n_season_compost", "n_season_lime", "n_season_fallow")
cor(d[,inputUse], use="complete.obs")

                 n_season_fert n_season_compost n_season_lime n_season_fallow
n_season_fert       1.00000000       0.35646933    0.17075604     -0.06439632
n_season_compost    0.35646933       1.00000000    0.03828723     -0.16923273
n_season_lime       0.17075604       0.03828723    1.00000000     -0.01116383
n_season_fallow    -0.06439632      -0.16923273   -0.01116383      1.00000000

Interpretation: The strongest correlation between the input use intensity variables is between seasons of fertilizer and compost use, ~0.35. While this is on the higher end it’s not necessarily cause for alarm.

inputUse <- paste(c("n_season_fert","n_season_compost", "n_season_lime", "n_season_fallow"), collapse= " + ")
list1 <- lapply(soilVars, function(x){
  mod <- lm(as.formula(paste("d[,x] ~", inputUse, "+ as.factor(cell)", sep="")), data=d)
  return(mod)
})
stargazer(list1, type="html", 
          title = "2016A Rwanda Soil Health Baseline - Naive Agronomic Practice Models",
          covariate.labels = c("Seasons of Fertilizer", "Seasons of Compost", 
                               "Seasons of Lime", "Seasons of Fallow"),
          dep.var.caption = "",
          dep.var.labels = "",
          column.labels = c(gsub("m3.","", soilVars)),
          notes = "Includes FE for cell",
          omit=c("cell"), out=paste("output", "rw_baseline_agprac.htm", sep="/"))

2016A Rwanda Soil Health Baseline - Naive Agronomic Practice Models

	Ca	Mg	pH	Total.N	Total.C	ExAl
	(1)	(2)	(3)	(4)	(5)	(6)
Seasons of Fertilizer	-5.326	-1.128	-0.001	-0.0004**	-0.003	-0.001
	(3.321)	(0.702)	(0.003)	(0.0002)	(0.003)	(0.004)
Seasons of Compost	1.727	0.363	0.004	0.0001	0.002	-0.006**
	(2.680)	(0.567)	(0.003)	(0.0001)	(0.003)	(0.003)
Seasons of Lime	41.659***	7.998***	0.024*	0.002**	0.022	0.009
	(14.473)	(3.061)	(0.014)	(0.001)	(0.015)	(0.016)
Seasons of Fallow	-19.069***	-3.517***	-0.028***	0.0005	0.009*	0.025***
	(5.283)	(1.117)	(0.005)	(0.0003)	(0.005)	(0.006)
Constant	548.378	68.041	4.684***	0.199***	3.360***	1.904***
	(387.623)	(81.977)	(0.370)	(0.021)	(0.392)	(0.420)
Observations	2,443	2,443	2,443	2,443	2,443	2,443
R2	0.555	0.585	0.612	0.502	0.456	0.570
Adjusted R2	0.514	0.547	0.577	0.456	0.407	0.531
Residual Std. Error (df = 2238)	387.475	81.946	0.370	0.021	0.392	0.420
F Statistic (df = 204; 2238)	13.684***	15.452***	17.322***	11.044***	9.202***	14.545***
Note:	p<0.1; p<0.05; p<0.01
	Includes FE for cell

Interpretation The naive model suggests that when we include site level fixed effects, duration of agronomic practices don’t have a big effect on soil health outcomes. However, some of the practice intensity variables are not well distributed. Let’s take a look at a log transformation. I’m adding one to the variables as to not end up with lots of Inf values.

Log transformations in theory are appropriate for variables that are right skewed (vavlues clustered to the left of the distribution) and see diminishing returns to increasing values. The shape of the data suggests a log transformation but it’s debateable whether the relationship is diminishing.

agPrac <- c(names(d[grep('n_season_', names(d))]))
for(i in 1:length(agPrac)){
  print(
    ggplot(d, aes(x=d[,agPrac[i]])) + geom_density() +
      labs(x = paste(agPrac[i], " No transform", sep = ""))
        )
}

# since these are all skewed, consider a log transform 
for(i in 1:length(agPrac)){
  print(
    ggplot(d, aes(x=log10(d[,agPrac[i]]+1))) + geom_density() + 
      labs(x = paste(agPrac[i], " Log transform", sep = ""))
        )
}

# look at other transfomations
for(i in 1:length(agPrac)){
  print(
    ggplot(d, aes(x=d[,agPrac[i]]^(1/3))) + geom_density() +
      labs(x = paste(agPrac[i], " cubic root transform", sep = ""))
        )
}

# visualize the outcomes as well to see if a transformation is warranted
for(i in 1:length(soilVars)){
  print(
    ggplot(d, aes(x=d[,soilVars[i]])) + geom_density() +
      labs(x = soilVars[i], title = soilVars[i])
        )
}

d$logFert <- log(d$n_season_fert+1)
d$logCompost <- log(d$n_season_compost+1)
d$logLime <- log(d$n_season_lime+1)
d$logFallow <- log(d$n_season_fallow+1)

Or look at BoxCox graph to empirically determine the right transformation. Log is assuming a diminishing return to an increasing X. That’s probably not the case with fertilizer. We’d actually expect an increasing return as values get larger. We use boxcox to see what the data suggest. We interpret it as follows:

• lambda = 2 -> square
• lambda = 1 -> no transformation
• lambda = 0.5 -> square root
• lambda = 0 -> log
• lamdba = -1 -> inverse

library(MASS)
for(i in 1:length(agPrac)){
boxcox(lm(pH ~ d[,agPrac[i]], data=d))
}

For pH at least it seems like a log transform is appropriate. We can run this for all other variables as well to see what we get as well.

Let’s look at the log results: See here and here for guidance on intepreting log transformed right hand side variables. See here for additional guidance on choosing a transformation.

How to interpret RHS log transform: For a linear multivariate OLS regression, we say “a one unit increase in X causes a (coefficient) change in Y.” For a linear-log regression where the X variable is log transformed, we say a L percent change in X leads to a (coefficient*L) change in Y.

Question: Do we want to run this analysis for only OAF farmers? If so, adjust the sample accordingly.

logVars <- paste(names(d[grep("log", names(d))]), collapse=" + ")
list2 <- lapply(soilVars, function(x){
  mod <- lm(as.formula(paste("d[,x] ~",  logVars, "+ as.factor(cell)", sep="")), data=d)
  return(mod)
})
list2b <- lapply(soilVars, function(x){
  mod <- lm(as.formula(paste("d[,x] ~",  logVars, sep="")), data=d)
  return(mod)
})

suppressWarnings(
  stargazer(list2, list2b, type="html", 
          title = "2016A Rwanda Soil Health Baseline - Log Agronomic Practice Models",
          covariate.labels = c("Seasons of Fertilizer (log)", "Seasons of Compost (log)", "Seasons of Lime (log)", "Seasons of Fallow (log)"),
          column.labels = c(rep(gsub("m3.","", soilVars),2)),
          dep.var.caption = "",
          dep.var.labels = c("",""),
          add.lines = list(c("Cell FE?", rep("Yes", 5), rep("No", 5))),
          notes = "Includes FE for cell",
          omit=c("cell"), out=paste("output", "rw_baseline_agprac_log.htm", sep="/"))
)

plm.log <- function(x, range){
  
  beta = round(summary(x)$coefficients[range,1],3)
  beta.pval = summary(x)$coefficients[range,4]
  beta.conv = ifelse(beta.pval < 0.01, "***", ifelse(
    beta.pval < 0.05, "**", ifelse(
      beta.pval < 0.1, "*", "")))
  #beta.pval = round(beta.pval, 3)
  outcome = paste(beta, beta.conv, sep = "")
  outcome = c(outcome, unique(round(summary(x)$coefficients[1,1],3)))
  res = data.frame(outcome, stringsAsFactors = F)
  
  return(res)  
}
rw.table15 <- do.call(cbind, lapply(list2, function(x){
  plm.log(x, 2:5)
}))
colnames(rw.table15) <- soilVars
rownames(rw.table15) <- c(names(d)[grep("log", names(d))], "constant")

4.2 Table 15 outcomes

rw.table15 <- rw.table15[,c("pH", "Total.C", "Total.N", "m3.Ca", "m3.Mg", "ExAl")]
write.csv(rw.table15, file=paste("output", "rwtable15_reg.csv", sep = "/"),
          row.names = T)

# a 10 percent increase in x leads to B1*.1 change in Y
logTrans <- do.call(cbind, lapply(list2, function(x){
  coeff = x$coefficients[2:5]
  tenPercent = round(coeff*.1, 5)
  names(tenPercent) <- paste("10% increase in ", gsub("log","", names(tenPercent)), " leads to:", sep="")
  return(tenPercent)
}))
colnames(logTrans) <- soilVars

print(kable(logTrans))

	m3.Ca	m3.Mg	pH	Total.N	Total.C	ExAl
10% increase in Fert leads to:	-2.31951	-0.52268	-0.00091	-0.00016	-0.00132	-0.00013
10% increase in Compost leads to:	1.00464	0.17075	0.00180	0.00007	0.00074	-0.00217
10% increase in Lime leads to:	12.28263	2.25366	0.00658	0.00053	0.00676	0.00040
10% increase in Fallow leads to:	-7.04768	-1.33558	-0.01025	0.00011	0.00222	0.00895

Interpretation: When we transform the variables to log, the data starts to tell a more coherent story, at least directionally. If we remove the district FE, the coefficients gain significance.

• Additional seasons of fertilizer use relates to a decrease in soil N and C.
• Seasons of compost and lime have the opposite, positive effect on soil N.
• Unsurprisingly, soil health metrics are poorer for soils on which more fertilizer has been used and lime has not been used.

4.2.1 Tenure with One Acre Fund

Let’s look first at a naive model of One Acre Fund tenure on soil health. Remember: these data are not longitudinal! These data are not longitudinal and reflect farmer selection into One Acre Fund. While these models will try to be both robust and parsimonious, we will inevitabily suffer omitted variable bias due to a lack of an instrument.

list3 <- lapply(soilVars, function(x){
  mod <- lm(as.formula(paste("d[,x] ~",  "total.seasons + as.factor(cell)", sep="")), data=d)
  return(mod)
})
stargazer(list3, type="html", 
          title = "2016A Rwanda Soil Health Baseline - Naive Tenure Models",
          covariate.labels = c("OAF Tenure"),
          dep.var.caption = "",
          dep.var.labels = "",
          column.labels = c(gsub("m3.","", soilVars)),
          notes = "Includes FE for cell",
          omit=c("cell"), out=paste("output", "rw_baseline_tenure.htm", sep="/"))

2016A Rwanda Soil Health Baseline - Naive Tenure Models

	Ca	Mg	pH	Total.N	Total.C	ExAl
	(1)	(2)	(3)	(4)	(5)	(6)
OAF Tenure	-0.817	-0.104	0.002	-0.0002	-0.001	-0.006*
	(2.856)	(0.603)	(0.003)	(0.0002)	(0.003)	(0.003)
Constant	556.919	69.700	4.693***	0.200***	3.370***	1.893***
	(389.276)	(82.243)	(0.373)	(0.021)	(0.392)	(0.422)
Observations	2,443	2,443	2,443	2,443	2,443	2,443
R2	0.550	0.581	0.606	0.500	0.455	0.565
Adjusted R2	0.510	0.544	0.570	0.455	0.406	0.527
Residual Std. Error (df = 2241)	389.234	82.235	0.373	0.021	0.392	0.422
F Statistic (df = 201; 2241)	13.648***	15.480***	17.114***	11.128***	9.300***	14.510***
Note:	p<0.1; p<0.05; p<0.01
	Includes FE for cell

Interpretation: The naive One Acre Fund tenure model suggest that across the board that additional years of 1AF practices have a negative effect on soil health parameters. Let’s combine 1AF tenure with the agronomic practices model above to build a more robust model:

Interpretation: Including agronomic practices and 1AF tenure in the same model dampens the magnitude, but not the significance, of 1AF tenure on soil health outcomes.

4.3 Agronomic practices - 15B cultivation practices

Thus far we have looked at aggregated historical plot level practices and their effect on soil health. We also asked farmers about their cultivation practices on their plot in the previous season, 15B. We have more precise information for fertilizer, compost and liming practices for the 15B season.

4.3.1 Graph the application/are variables

Conclusion - Take 1: The application rate per are variables are weird. I think it’s because of the field dimensions. I’m going to go back to the field dimensions and check this.

Let’s look at the dimensions of the fields that have large application rates

4.4 Tenure Scatter plots

Generate tenure vs. pH and tenure vs. fertilizer seasons scatter plots:

4.5 Table 16 - Kgs of soil ammendments - fertilizer (kg)

4.6 Table 16 - all fertilizer used

4.7 Table 16 - Kgs of soil ammendments - compost (kg)

The compost regression results are suprising. Look at Carbon and compost scatter plot. If we include 0 compost, we see no relationship. If we exclude 0 compost, we

There are too many zero or near zero values. We can conclude that perhaps we’re not asking this question well. Even when we remove non-zero values, there’s little discernable relationship.

4.8 Table 16 - Kgs of soil ammendments - lime (kg)

4.9 Table 16 - types of fertilizer - DAP (kg)

4.10 Table 16 - types of fertilizer - NPK (kg)

4.11 Table 16 - types of fertilizer - Urea (kg)

stargazer(list5, type="html", 
          title = "2016A Rwanda Soil Health Baseline - 15B practices",
          covariate.labels = c("Fertilizer Rate (log)"),
          dep.var.caption = "",
          dep.var.labels = "",
          column.labels = c(gsub("m3.","", soilVars)),
          notes = "Includes FE for cell",
          omit=c("cell"), out=paste(od, "rw_baseline_15b_ag.htm", sep="/"))

4.12 Farmer fertility perception

Let’s look at farmer perceived fertility as a predictor of soil health. We’ll set ‘same fertility’ as the reference category.

Interpretation: Farmers understand their fields well. Their categorization of which field are more and less fertile corresponds to our quantified measures of soil health. The only features farmers don’t seem to get correct are nitrogen and carbon. The nitrogen and carbon levels are indistinguishable in ‘low fertility’ fields relative to the fields deemed to be the ‘same fertility.’ Reminder: We need to remember that farmers are only evaluting one of their fields thus we are not able to account for the quality of the farmer in assessing his/her fields.

4.12.1 Deeper dive into farmer perception of fertility

• Look at wet chem values for all soil features by farmer perception
• Run PCA on all soil attributes to see which principle components predict soil fertility
• Consider adding additional control variables into model (like Ca as control for pH perception model)

4.12.2 For Eric, CEC summary for new Ca and Mg thresholds

Interpretation: Our sample size decreases considerably when we look only at wet chemistry results. Thus, we do not see the same signficant relationships we saw between farmer perceived fertility and soil characteristics we saw when looking at the predicted values.

4.13 Principal Component Analysis (predicted values)

The first principal component is composed primarily of soil pH related variables, Ca, Mg, and pH. The second capture N and C. These groupings (pH grouping and C and N) are not all that surprising given that our predicted soil variable set is fairly limited.

If the variables have the same sign that indicates that they are positively correlated in the principal component. In the first principal component, we see Ca, Mg and pH loading in the same direction. In the second principal component, we see N and C moving in the same direction. Ca and Mg are also associated in the same direction but to a lesser extent.

When we layer farmer perceived soil fertility on top of the principal component scatter plot, no clear pattern emergres. Fields with the same fertilitiy, less and more fertility are indistinguishable by the principal components. Thus, there is not a clear profile for farmer identified healthier soil or weaker soil based on principal components.

4.13.1 Geographic Descriptors - PCA for soil profiles

This section will build on the principal component work above and look at improving understanding of local context to inform local adaptation. Here we will also construct soil profiles to simplify the scaling of promising products and practices to targeted locations. We don’t have a full suite of predictors so we can’t look at a comprehensive soil profile.

When we color the figure by district, we start to see a pattern emerge. However, there are too many points to clearly detect where all districts fall. Let’s instead look at the data by AEZ.

Coloring the points by AEZ, we see a much clearer trend. The east is to the right of our graphic, then the central plateau followed by lake Kive and then Congo Nile in green on the left. What does this mean in terms of actual soil features? Let’s look at the figure but with the variable associations layered on top. See here for documentation

It appears that the right side of the graph, the eastern AEZ, is associated with pH, Mg and Ca. This indicates that as the first principal component increases, so does the level of pH, Ca and Mg. Conversely, total N and C increase as the second principal component decreases. In terms of the AEZ figure above, this suggest that the Congo Nile AEZ has higher levels of N and C. Let’s test these hypotheses with a simple summary table:

Confirmed! This likely also suggests that Congo Nile is at a higher altitude than the surrounding areas and that eastern Rwanda has relatively less weathered soils compared to western.

Consequence for trial placement: The Rwanda program is already blocking trials by AEZ but these data confirm that AEZ reflects meaningfu soil variation and thus captures key growing conditions for Rwandan farmers. Blocking trials by AEZ in Rwanda will enable us to evaluate trial hypotheses in more neutral pH ranges and higher N and C conditions.

4.13.2 Principal Component Analysis (wet chemistry values)

4.14 District and cell level summaries of soil and managment practices

4.15 Female farmers farming poorer soils?

4.16 Export the data for Eric

5 Propensity Score Matching

We need to do a more rigorous job of accounting for differences between Tubura farmers and identified control farmers. Execute propensity score matching (PSM) to identify control farmers that overlap with Tubura farmers with regard to their likelihood of being a Tubura farmer.

Interpretation We have some overlap but it’s clear that Tubura farmers occupy a different range than the identified control farmers. Let’s continue with the PSM matching process but restrict ourselves to Tubura and control farmers that meet a certain PSM matching radius.

5.0.1 Notes on PSM process:

We have to indicate a variable for matching. I’m choosing pH as we know it to be an issue in most of our operating areas and addressing soil acidity has numerous residual benefits to soil health. I’ll want to do this for all soil outcomes, however.

• As a next step, I’d like to identify the matched subset and then regress Tubura tenure on soil outcomes. Does that make sense?
• I chose a caliper of 0.25. I should review this process with Maya to make certain I’m following best practice. I sort of just pulled that sd figure out of the sky.
• An initial check of the match quality using a common model indicates that the matches are poor.

Now check the naive model approach for PSM balance.

Interpretation: We want to see standard mean differences less than the absolute value of 0.25 (or 0.1 if we’re being conservative) and variance ratios close to 1 but certainly between 0.5 and 2.

According to the CRAN summary, sdiff is the standardized difference between the treatment and control units multiplied by 100. If I divide by 100, the values come much closer to reasonable value.

5.0.2 Let the PSM models vary

The common model approach doesn’t seem to be working for any of the variables. I’m going to rework the modeling approach so we can fit different models for each outcome upon which we’re trying to match.

The models can now vary by outcome. Let’s see if we can improve our results.

Interpretation If I divide the standardized mean differences by 100, we meet the balance criteria of the standardized mean difference being close to 0 and the variance being close to 1. Let’s print out the model results to see how Tubura and control farmers compare on key soil meterics. These results should supercede the naive balance tables presented above.

We achieve acceptable balance for the soil attributes but we don’t for seasons of fertilizer use. This is isn’t entirely unexpected given that Tubura’s primary service is providing fertilizer inputs and training.

I’m comparing season 1 clients to clients with more tenure in the PSM matching but I also need this variable for the threshold calculations

5.1 Calculate farmers under soil health thresholds - table 1

New clients vs. returning clients summaries:

New Clients vs. returning clients thresholds

5.2 Table 1:3 coefficients

Interpretation: Propensity score matching gives us a comparable treatment and control group. The table above shows that after matching on those characteristics, there are effectively no differences between One Acre Fund farmer and Tubura farmers on soil attributes. The unadjusted p-values show 1AF farmers to have slow levels of soil nitrogen but this finding disappears if we account for running multiple matching models.

When we expand the outcome variable set to include practice variables, we first no longer get a good propensity score match for all variables.

• Tubura participation so strongly overlaps with fertilizer use that the comparison group is no longer comparable. The resulting model shows a significant difference in seasons of fertilizer use.
• Tubura farmers have used compost for 0.3 more seasons than non-Tubura farmers. We don’t observe a difference in N or C levels in the soil attributes however.
• We don’t see a difference between Tubura and non-Tubura farmers on seasons of legumes being the primary crop in the study plot.

5.3 Break out first season clients from more tenured - for appendix?

This requires us to re-run the PSM matching process, confirm that our models are sound and then output the results ala table 2 again.

5.4 Study group balance post-match

See here for some guidance on hwo to use weights to reconstruct the group balance following the matches.

See here for weighted t.test documentation

5.5 Models with PSM matched farmers

Per Robert’s suggestion, now that we’ve matched Tubura and non-Tubura farmers, let’s assess the severity of Tubura tenure on key soil health outcomes.

suppressWarnings(
stargazer(tenureTab.add, type="html", 
          title = "2016A Rwanda Soil Health Baseline - PSM Tenure",
          covariate.labels = "One Acre Fund Tenure", 
          dep.var.labels = "",
          #column.labels = modNames,
          #column.labels = c(gsub("m3.", "", as.vector(namesInput))),
          notes = "Includes FE for cell",
          omit=c("cell"), out=paste(od, "rw_baseline_matched_tenure.htm", sep="/"))
)

Interpretation: Using a PSM matched sample, the models above assess the effects of additional years of farming with Tubura. Numerous control farmers have also been Tubura farmers in previous seasons. Thus, I’m keeping the model simple instead of adding a client*tenure interaction. We can easily test that as well though.

6 Soil Clusters

6.1 Heirarchical Clustering

I draw from Robert’s example in the BGMS soil response notebook.

Scale soil variables to remove units as an issue.

Note: I’m honestly not certain what we’re looking for in these outputs. Look for a balanced tree? Ward looks pretty balanced.

6.2 K-means clustering

K-means clustering is a typical unsupervised learning algorithm. We start by identifying the number of clusters we should tell the formula to look for.

We’re looking for the bend in the graph. It seems to come at 3 but we can try both 3 and 4 and see how they match with our understanding of AEZ and environmental conditions.

These clusters are roughly equally sized. Robert used another method, partitioning around Mediods, to better account for outliers that were causing the algorithm to create a small additional cluster.

Give the clusters brief descriptions

6.3 Mappaing the location of the clusters

6.4 Soil clusters by altitude and AEZ

6.4.1 Soil clusters by altitude and rainfall

Try this when we have aWhere data

6.5 Compare clusters to AEZ

We’re not seeing close alignment between clusters and AEZ. This is probably because AEZ are a mix of altitude, rainfall, temperature and soil. When we have weather data we can try to cluster again and see if we’re more closely aligned.

Alternatively, this could be initial evidence that the AEZ designations are not capturing and separating variation very well. I think we need to start with more data though.

7 Mapping

7.1 Soil health maps for soil health study areas

Generate district density summaries for Nathaniel to pair with the map

7.2 Static map of sampling density

7.3 Plot simple summary of soil characteristics

Consider revising these maps to a smaller greographic unit. Add the name of the location for uninitiated users.

# This is a small experiment to combine raster (spdf) and leaflet and be able to access the data in the raster interactively.
#mapLayer <- sp::merge(rw, ss.soil, by.x="NAME_3", by.y="Group.1")

# fill = T, fillOpacity = 0.7, fillColor = d.fill, 
#         stroke = T, color = "white", weight = 2, dashArray = 3, 
#         opacity = 0.5, popup = county.tt(d)

# leaflet(mapLayer) %>% addTiles() %>%
#   setView(lng=rwanda$longitude, lat=rwanda$latitude, zoom=8) %>%
#   addPolygons(
#     stroke = TRUE, opacity=0.2, smoothFactor = 0.5, 
#     fillColor=mapLayer$pH, fillOpacity = 0.5)

7.4 Interpolations of soil health values

Interpolate soil health values for full operating area using soil health study values. We want to eventually add all Rwandan soil values into a single dataset to update and hone these values. See here for more guidanace

note:

• Layer on the Tubura sites to the interpolated map.
• Use leaflet so you can zoom in and out more easily.
• Make it a raster layer that you can click on understand the values at different locations and the name of the location. I’ll need to make this a shiny app to have that functionality

The code below will run 5 K-fold cross validation to compare interpolation models. The output will be fed into the interpolate leaflet code below.

Check that I’m handling the projection correctly with Robert

Now loop over the variables of interest where x is the soilVar variable.

7.5 Interpolated soil maps

Print out the interpolated values for inclusion in the report.

–end